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ABSTRACT 


A new three-dimensional numerical cloud model has been developed for the 
general purpose of studying convective phenomena. The model utilizes a time 
splitting integration procedure in the numerical solution of the compressible 
nonhydrostatic primitive equations. Turbulence closure is achieved by a 
conventional first-order diagnostic approximation. Open lateral boundaries 
are incorporated which minimize wave reflection and which do not induce 
domain-wide mass trends. Microphysical processes are governed by prognostic 
equations for potential temperature water vapor, cloud droplets, ice crystals, 
rain, snow, and hail. Microphysical interactions are computed by numerous 
Orville-type parameterizations. A diagnostic surface boundary layer is 
parameterized assuming Monin— Obukhov similarity theory* The governing 
equation set is approximated on a staggered three-dimensional grid with 
quadratic-conservative central space differencing. Time differencing is 
approximated by the second-order Adams— Bashf or th method. The vertical grid 
spacing may be either linear or stretched. The model domain may translate 
along with a convective cell, even at variable speeds. In storm splitting 
cases, the domain translates with the convective cell having cyclonic rotation 
and allows the other cell(s) to pass through the lateral boundary without 
detrimental consequences . 

Potential applications of the model range from the simulation of shallow 
cumulus to supercell cumulonimbus, including such convective phenomena as 
downbursts, tornadoes, gust fronts, and hailstorms. 
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1 . INTRODUCTION 


Three-dimensional, convective cloud models have now advanced to a stage 
where they can be directly compared to observed data fields. This has been 
made possible by the current evolution of high-speed and large in-core memory 
vector computers. The earliest three-dimensional cloud models were developed 
by Steiner (1973), Miller and Pearce (1974), Pastushkov (1975), Schlesinger 
(1975, 1978), Lipps (1977), Klerap and Wilhelmson (1978a), Cotton and Tripoli 
(1978), and Clark (1979). These pioneering models were limited by computer 
restraints, and were run with relatively crude grids and simple 
microphysics. Nevertheless, they were able to produce much information on the 
dynamics of buoyant convection within vertically-sheared environments. 
Refinements in these models have continued to progress (e.g., Yau, 1980; Cho 
and Clark, 1981; Wilhelmson and Chen, 1982; Tripoli and Cotton, 1982; Yau and 
Michaud, 1982; Schlesinger, 1984a, 1984b; Smolarkiewicz and Clark, 1985); but 
only a few, so far, have attempted to verify their model simulations with 
detailed observed data sets. Of the 3-D models, only Cotton et al. (1982) 
have included parameter izations of ice-phase microphysics. They found that 
its inclusion moderately affected the dynamics of the simulated clouds. The 
significance of including the ice phase has also been shown in other studies 
with one- and two-dimensional models. For example, Ogura and Takahashi (1971) 
have found that the exclusion of the ice-phase resulted in a considerable 
change in the evolution of the downdraft. Willoughby et al. (1984) and Lord 
et al. (1984) using a 2-D axisymmetric model, have found that inclusion of 
ice-phase microphysics resulted in dramatic differences in a hurricane 
simulation; one important finding was that the locations of mesoscale 
downdrafts were controlled by falling ice particles. Ice-phase microphysics 
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cannot be casually neglected from model developments. It may have an 
important impact, especially with regard to simulations of downburst phenomena 
and deep tropospheric convection. 

The purpose of this report is to present the development of a new three- 
dimensional numerical model. The model, the Terminal Area Simulation System 
(TASS), has a meteorological framework and is formulated for the general 
purpose of studying the physical-dynamical character of convective clouds and 
storms. The TASS model is capable of realistic simulations of convective 
clouds ranging from nonprecipitating cumulus to intense, long-lasting, 
supercell hailstorms. Its application, however, is not limited to convective 
clouds; the model may be applied to many other microscale and meso-gamma scale 
phenomena. In fact, considerable care has been taken in the formulation, so 
that the model is capable of valid simulations of tornadic and severe 
downburst phenomena. One major use of the model, so far, has been in the 
study of downburst-r elated wind shear and its impact on aviation safety (e.g., 
Chuang et al., 1984; Proctor, 1985a, 1985b). The model is currently being 
used to examine the three-dimensional structure of downbursts and to provide 
realistic data for real-time flight simulations. 

A brief description of the TASS model is as follows. The model utilizes 
a nonhydrostatic, compressible and unsteady set of governing equations which 
are solved on a three-dimensional staggered grid. The model divides water 
into six bulk categories each governed by a prognostic equation. The six 
categories are: 1) water vapor, 2) ice crystals, 3) cloud droplets, 4) rain, 

5) snow, and 6) hail/graupel. The former three categories represent 
nonprecipitating forms of water, while the latter three represent 
precipitating forms of water. The hail/graupel category may consist of either 
hail or graupel. Note that all three phases of water (i.e., vapor, liquid, 
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and solid) are included. Parameterization of the numerous microphysical 
interactions (that result in exchanges of water between the six categories) 
are similar to those given in Lin et al. (1983), and Rutledge and Hobbs 
(1983). As for treating turbulence mixing, the TASS model adopts the subgrid 
closure approach (e.g., Deardorff, 1970; 1972; 1973). That is, scales of 
turbulence larger than the assumed grid size are simulated explicitly within 
the flow field; while scales of turbulence less than the grid size are 
parameterized from a closure approximation. The subgrid closure model 
currently in use is a conventional, first-order, diagnostic approximation. 

TASS also incorporates surface stresses which are dependent upon 
stratification, ground roughness, and local winds. Numerical stability and 
conservation in the solution of the governing equations relies on an 
appropriate choice of numerics and boundary conditions. The TASS model uses 
quadratic— conservative space differencing and incorporates a modified Orlanski 
radiation boundary scheme. Application of the radiation boundary condition to 
the open lateral boundaries allows the outward propagation of waves with 
minimal reflection. Also, the procedure for applying the radiation boundary 
conditions is free of domain-wide mass trends. Other features of TASS are 
1) the option of a vertical grid-size stretching, 2) movable mesh with time 
varying translation speed, 3) a numerical filter and sponge applied below the 
top boundary, and 4) specification of an initial environment from a sounding 
that is either observed or predicted from a regional model simulation. Output 
from TASS includes three-dimensional fields of wind velocity, rain, snow, 
hail, cloud water (cloud droplets and ice crystals), radar reflectivity, 
temperature, and pressure (see Fig. 1). 

Details of the model formulation are found in Chapters 2-6. In Chapter 2 
the basic model assumptions are listed. In Chapter 3 the model framework. 
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governing equations, boundary conditions, and turbulence closure are 
discussed. The cloud microphysics, including the development of the 
microphysical parameterizations, are described in detail in Chapter 4. The 
initial and reference conditions are discussed in Chapter 5. This section 
includes the formulation for the initial perturbation fields which are 
necessary in order to trigger convective development. In Chapter 6 the 
numerical procedure is described. This section includes the formulation for 
the variable-speed, grid-translation algorithm, as well as the details of the 
finite-difference equations, grid, and numerical stability criteria. Also 
included in Chapter 6 are some important numerical details in the computation 
of the cloud microphysics and a brief discussion of the model code. 

Several test simulations with the TASS model are described in Chapter 7 . 
These test cases assume simplified atmospheric conditions, and are useful in 
demonstrating the validity of the model coding and formulation. A more severe 
test of the model performance is discussed in a second report, in which TASS 
simulated results are compared and evaluated against observed data sets. 
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Fig. 1. Schematic of Data Flow 







2. BASIC MODEL ASSUMPTIONS 


The TASS model has been developed for the general purpose of studying 
convective phenomena with time scales of several hours or less. The primary 
model assumptions are: 

1) an equation set which is valid for subsonic and high-Reynold' s number 
turbulent flow; 

2) Reynolds' averaging of equation set is roughly equal to grid size; 

3) thermal radiation is neglected; 

4) the first grid point above the ground lies within the surface stress 
layer ; 

5) the ground is flat with a homogeneous surface roughness; 

6) only a passive interaction with the large-scale environment 
disturbances can propagate out of the limited model domain, but not 
into the domain; 

7) the initial environment is horizontally homogeneous and in steady 
balance — convection is initiated by adding a velocity and/or 
temperature perturbation; 

8) super saturation with respect to liquid water is not allowed — 
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condensation occurs at a rate which maintains saturation; 


9) subgrid-scale condensation is neglected; 

10) hydrometeors are classified into five bulk categories and 
microphysical interactions are parameterized; 

11) rain, snow, and hail/graupel assume inverse-exponential size 
distributions; 

12) cloud ice crystals have a monodispersive size distribution; 

13) falling hydrometeors instantaneously achieve their terminal velocity 
and have no horizontal slip relative to air motion; and 

14) electrical effects (e.g., drop charging) are ignored. 
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3. DYNAMIC MODEL 


Model Framework 

The model framework assigns a reference environment which is a function 
of only the height coordinate and is in hydrostatic balance. The dependent 
variables may be expanded in terms of the reference environment as (symbols 
are listed in the Appendix A) 

u(x,y,z,t) = U Q (z) + u'(x,y,z,t), 

v(x,y,z , t) = V Q (z) + v'(x,y,z,t), 

P(x,y,z,t) = P Q (z) + P(x,y,z,t), 

9(x,y,z,t) = 9 q (z) + 9'(x,y,z,t), 

p(x,y,z,t) = p Q (z) + p'(x,y,z,t), 

Q v (x,y,z,t) = Q v0 (z) + Q^(x,y,z,t), 

where x, y, and z are respectively the Cartesian coordinate in the west to 
east, south to north, and vertical direction; t is the time coordinate; u and 
v are respectively the west to east and south to north velocity component; P 
is pressure, 0 is potential temperature; p is the air density, Q v is the 
mixing ratio for vapor; U Q , V Q , P Q , p Q , 0 q , Q vq are the reference components, 
and u', v', p, p', 9' and are the deviations from their respective 
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reference quantities. The values of the reference environment may be taken 
from an actual hydrostatically-balanced environmental sounding. 

The reference environment also represents the initial model 
environment. Details of the initial and reference environment formulations 
are discussed in section 5. 


Governing Equations 

The model incorporates the unsteady primitive equations in nonhydrostatic 
and compressible form. As in Cotton and Tripoli (1978), dimensional pressure 
and potential temperature (along with moisture substances) are chosen as the 
prognostic thermodynamic variables. Closure of the equation set is obtained 
by diagnosing density and temperature. 

In deriving the following equations the hydrometeors are assumed to 
instantaneously achieve their respective terminal velocities; and thus, the 
total mass per unit volume of atmosphere is equal to the sum of the masses 
(per unit volume of atmosphere) of dry air, water vapor, cloud droplets, cloud 
ice crystals, rain, snow, and hail (see Appendix B). This assumption leads to 
an equation of state having the form 

p = If (1 “ °- 61 Q V + V’ (1) 

where Qj = Q CD + + Qr + QgN + Qh* Here, R is the gas constant for dry 

air; T is the air temperature; and Q T is the total of the water substance 
which is a sum of the mixing ratios of cloud droplets, Q CD , ice crystals, Q IC , 
rain, Q R , snow, Q SN , and hail, Q H . 

The governing equations consist of eleven prognostic equations. They are 
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expressed in a form that is consistent with the fully-elastic mass-continuity 
equation (i.e., dp/dt = — p V • V) as follows. 


Three equations for momentum 

The equations for momentum are expressed in Cartesian tensor notation as 


5u i + H_ dp = 

5t p dx. 

O 1 


5u,u, 5u 


2 Q u ' e . .. + — -r — ^ + (rr~) 

j k ijk p Q dXj dt 


( 2 ) 


The Einstein summation convention is used for vector quantities; u^ is the ith 
velocity component (i, j,k index from 1 to 3, u^ = u, u 2 = v * u 3 = w > x l = x > 
x 2 = y» and x 3 = z), is the alternating unit tensor, is the jth 

component of the Earth's angular velocity, g is the gravitational 
acceleration, and 6 is the Kronecker delta. Eq. (2) is in nonBoussinesq form 
as derived in Proctor (1982). The advection terms, represented by the first 
two terms on the right-hand-side of Eq. (2), are expressed in this expanded 
form in order that quadratic-conservative numerical formulations can be 
applied^ . 

Coriolis effects are retained even though Klemp and Wilhelmson (1978b) 
found that their inclusion produced only minor changes in cumulonimbus 
simulations. An approximate formulation of the coriolis effect follows that 
of Tripoli and Cotton (1982). They assume that the Initial fields are in 
geostrophic balance, and they neglect the horizontal variations of the initial 


1 The advection term can also be expressed as [3u u p /9x + u^u p Q /3x ]/p q 

(e.g. Tripoli and Cotton, 1980). This form, as well as the form expressed in 
Eq. (2) were tested in an axisymmetric simulation of a firestorm; even under 
these severe conditions results from both of the formulations were nearly 
identical. 
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pressure and temperature fields that are due to coriolis effects. These 
assumptions lead to the formulation in Eq. (2) in which the coriolis 
acceleration only affects the perturbation velocity. 

The next-to-last term in Eq. (2) is due to Reynolds' averaging. The 
Reynolds' stresses, T ±y are a result of subgrid-scale fluctuations of 
velocity. Details of the subgrid formulation are discussed later in this 
section. 

The last term in Eq. (2) is an external forcing term which is added in 
order to maintain a steady initial state. Details of this procedure will be 
discussed in section 5. 

The density ratio term H represents the ratio of the reference density of 
the environment to the local density. It may be diagnosed (see Proctor, 1982) 
as 


H = p o /p = (9/e o ) (P o /P) 1/r > [1 + 0.61 (Q v -Q vo ) - Q t ], (3a) 

where ri = C^/C^ the ratio of specific heats of air at constant pressure 
and constant volume. The exponentiated term in Eq. (3a) can be expanded, 
resulting in an expression which is more computationally efficient, yet still 
valid for most atmospheric problems (where p « P Q ). This alternate 
expression, which is used in almost all of the model experiments, is 

H = t 0/0 o " p/tiP o ] [1 -° + °* 61 ( V Q vo ) " Q t!* (3b) 


11 



Prognostic equation for pressure deviation 


The prognostic equation for pressure deviation is 


* 


j 


3u,p 5u 

Sip + p + P o gU 3 6 j3 


tip de . dq T dQ v 

+ r si' " p [ dr ' - 61 sr 


i <i - 


•61Q v + Q t ) 


+ (4) 

P o 

where S(P) is due to subgr id-scale fluctuations- The derivation of this 
equation can be found in Appendix B. 


Thermodynamic equation 

The thermodynamic energy equation is written for potential temperature 
which is conserved in dry~adiabatic processes. The prognostic equation for 
potential temperature is 


ae-_ &U J 9 + e !!i+ 

dt e>Xj ax^ p o axj 


* 



p 


where L v , Lf , and L g are respectively the latent heats of vaporization, 
fusion, and sublimation; S v , S f and S s represent the rate of phase conversions 
of water in units of mass of water per unit of time per unit mass of dry 
air. This basic formulation of the thermodynamic equation was derived by Das 
(1969) and has been found to have reasonable accuracy (e.g., Wilhelmson, 

1977). 
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Six equations for water subs tance 


A coupled set of six prognostic continuity equations govern the 
distribution of water vapor, and liquid phase and solid phase water 
substances. Each variable for water substance is expressed in terms of a 
mixing ratio (mass of water substance per mass of dry air) which is conserved 
in the absence of turbulence mixing, phase changes, and microphysical 
interactions. 

The prognostic equations for water vapor, nonprecipitating cloud 

droplets, nonprecipitating cloud ice crystals, rain, snow, and hail are, 
respectively, 


A 

at 

9Q 

at 

aQ 


au q au 

- r~ J + Q — J- + 

ax , v ax, p Sx 


CD 


at 

aQ 

at 

5Q, 

at 

aQ 


ic 


R 


SN 


j 

du j Q cp 

ax. 


au.Q. 


j "O ““j 
5u 


aS(Q ) 9Q * 

+ S + (— ^-) 

wfl •n 'A** ' > 


vap at 


+ Q 


du j A 1 95 . 


CD dx j p o 9x j 


Ac 5u i 1 as(Q ic> 

HP ^ • 


CD’ 


ax^ ' 1 C PXj P Q bXj ■ “ic* 

au. , aQ w P 


au j Q R 


bx 


+ q t 


j + 1 “' t R"R K o s L 1 


6„. + 


as(Q D ) 


j R ax j p o ax j ’ Po ax j 


+ S R‘ 


a u ,Q 


j Q SN . ft A . 1 5Q SN W S p o . , aS «W 

j v sn a Xj + p o a^ 6 3 j + + 


H 


at 


Su j q n +0 Su J + 1 8 q hVo . , 8S <V 

T Hr, — + X O. . + 


3x 


j 


h ax, p ax, 

j o j 


3j ax 


j 


+ V 


( 6 ) 

( 7 ) 

( 8 ) 
( 9 ) 


(10) 


( 11 ) 


The last term in Eqs. (6) - (11) is a source term resulting from microphysical 
interactions between each of the bulk categories. In the absence of 
evaporation from the ground, continuity of water substances requires: 
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0 . 


s 

vap 


+ S 


CD 


+ S IC + S R + S SN + S H = 


Since Eqs. (9) - (11) govern precipitating categories of water substance, 
a fall-out or slip term is included. In these terms the bulk terminal 
velocities for rain, snow and hail are represented by W R , W g , and W R , 
respectively. The terminal velocities are by definition positive and are 
directed vertically downward. Details of their formulation and the modeled 
microphysics will be discussed in section 4. 


Treatment of Subgrid Processes 


If the dependent variables in the governing equations are treated as 
averages over the grid volumes, the equations themselves should be treated as 
similarly averaged, giving rise to residual terms. In the presence of 
turbulence motion these terms represent subgrid Reynolds stresses in the 
momentum equation and subgrid eddy transport in the remaining prognostic 
equations. This approach to turbulence closure allows the resolvable eddies 
to be modeled explicitly, while the influence of the subgrid eddies (the 
effects of eddies approximately equal to or less than the grid size) are 
"parameterized" (e.g., Deardorff, 1970; 1972; 1973; Sommaria, 1976). An 
overview of subgrid scale modeling can be found in Herring (1979). 


Subgrid Reynolds stresses 

The subgrid Reynolds stress tensor according to lst-order closure theory 
(e.g., Clark, 1979) is 

T ij ' % *M V 
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in which the deformation tensor Djj is defined 


9 u i 9 u . 9u 

) = — + 1 _ £ !S r 

ij 9x ax 3 a x 0 i j ’ 


The subgrid eddy-viscosity for momentum is 


= (a A) 2 IdefI (1-R )°' 5 , 


( 12 ) 


where a is an empirical constant, A is the subgrid-turbulence length scale, 
IdefI is the absolute magnitude of the local rate of deformation, and R f is 
the Richardson flux number. The averaging length scale is related to the grid 


size as 


A = (2Ax • 2Ay • 2Az) 


1/3 


(13) 


Studies by Clark et al. (1977,1979) and Love and Leslie (1977) suggest that 
the averaging length should be related to twice the grid length rather than 
the grid length as in Deardorff (1978). The magnitude of the local rate of 
deformation is 


I™' - ^ » 2 


0.5 


)] 


where the divergence is defined as 


9u k 

<b = ? — - 

k 9x k 


(14) 


The local Richardson flux number is related to the local Richardson number 
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as 


for 


- 100 <_ R f 1 0.99, 


R f = (K T /K M )R i 

and the local Richardson number is defined (e.g., Duran and Klemp, 1983) as 


r 


R. = 

1 


iDEFl 


- h % + q sn + V 

30 

- — - + 0.61 7— (Q + Q rn + Q Tr - % n ) 
T 0 z 3z v CD IC "vo 


tor q C D + Q ic^°' 


3z ^CD + Q IC + q v^ 

- h (Q R + Q SN + Q h> 


for 0 C D + q ic > ° ; 


where the equivalent potential temperature is defined 


0 e " 9 exp [(L g Q v + L f Q IC )/C p T]. 

An arbitrary upper bound of 0.99 is enforced on R f in order to guarantee a 
minimal amount of diffusion. An arbitrarily lower bound of -100 is also 
assumed . 

In the above formulation, the subgrid eddy mixing is affected by the 
local shears through the deformation term, and is modified by the 
stratification through the Richardson flux number. When the stratification is 
neutral (R f = 0), Eq. (12) reduces to the Smagorinsky turbulence model which 
has been applied with great practical success in planetary boundary layer 
studies (e.g., Deardorff, 1970; 1972). Clark et al. (1977, 1979) has 
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examined several closure assumptions in simulations of homogeneous isotropic 
turbulence and has found that the Smagorinsky model adequately accounted for 
the transfer turbulence energy to the subgrid scales. They also have found 
that the optimal dissipation of grid-scale turbulence energy occurred when 
a - 0.186. This value falls within the range of other experimental and 
theoretical values, and should remain invariant of grid size (at least as long 
as A is within the inertial subrange). The impact of unequal grid sizes on 
this turbulence closure scheme is unknown. 

Subgrid eddy transport 

Also from lst-order closure theory, the subgrid covariances in Eqs. (5) - 
(11) are approximated as 

s ( q ) = dq/9Xj ( 15 ) 

where q = 0, Q v , Q^, Q IC , Q R , Q SN , or Q R . The subgrid eddy-mixing 
coefficient for heat, K<j>, is assumed for all scalar variables (except 

pressure). Its value is taken from theoretical considerations by Deardorff 
(1973) to be 

Kr = 2.5 f or z > Az . 

Similar to Eq. (15) the subgrid transport for pressure deviation is 
assumed as 


S(p) = p o k m a P /a Xj . 
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Boundary Conditions 


Surface boundary 

The choice of surface boundary conditions can have an important influence 
in convective cloud simulations. For instance, Schlesinger (1982) found that 
a significant impact on simulated storm development resulted, when he changed 
his model surface winds from the customary free-slip to semi -slip. The change 
to the semi— slip boundary condition was made so as to take into effect the 
retarding effect of the earth's surface. Schlesinger found that the 
orientation of the flanking line and the dominance of the right moving 
convective cell (after storm splitting) were strongly affected by the choice 
of surface boundary conditions. The ground boundary layer also was found to 
play a crucial role in the axisymmetric simulations of a tornado by Proctor 
(1982). His model simulations demonstrated that surface convergence induced 
by a parent vortex, and the subsequent release of latent heat of condensation, 
may lead to the formation of a tornado. His simulations also demonstrated 
that frictional convergence was responsible for the extreme upward velocities 
in a tornado at a few tens of meters above the ground. These studies and 
others have demonstrated that friction at the earth’s surface can exert some 
influence on convective systems. Friction at the earth's surface may alter 
the depth, speed, and orientation of storm-produced surface outflows, and may 
also influence the propagation and intensity of convective storms. 

In the absence of detailed surface-layer data below convective storms, it 
is difficult to formulate and test surface boundary conditions for cloud 
models. In this model as in Sommeria (1976), a constant flux or stress layer 
is assumed to extend from the ground to the first grid point above the surface 
(at height h = Az/2). Conditions within this layer are parameterized using the 
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nonditnensional shear and temperature gradient functions, as deduced from field 
observations by Businger et al. (1971). The model surface formulation 

described below is completely diagnostic; it does not include a moisture or 
temperature budget for the ground. 

At z = h: 

< t > °' 5 

(kh)2 l DEF l C 1 -^V . 


and 




where k is von Karman's constant (0.4), and and ^ are respectively the 
nondimens ional wind shear and temperature gradient. Within the surface layer 
(0 < z < h), relationships for u, v, 9, and 1^ are based on Monin-Obukhov 
similarity theory and are derived in Appendix C. The mean velocity gradients 
in the surface layer are 

" u < h > /h » and <§J> - v(h)/h. 

The mean temperature gradient in the surface layer is 

V 1 + * 61 Q vo > [u(h)2+ v ( h ) 2 ] G 

^ “ 3 ~2 

8h V 

where G R and G M are universal functions that can be deduced from field 
observations and L is the Monin-Obukhov length. The mean eddy-viscosity for 
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momentum within the surface layer is 


<v ■ v h) vv 


The values for 4> M and $ R are determined 
(1971), i.e., 

unstable (h/L) < 0 

cf> = [1 - 15 (h/L )] -0,25 
M 

<j> = 0.74 [1-9 (h/L )] -0 * 5 

H 


from measurements of Businger et al. 


stable (h/L) > 0 


=1 + 4.7 (h/L) 

M 

4> = 0.74 + 4.7 (h/L). 

H 


Likewise, the universal functions are 


G m = An (h/z Q ) - 

G r = 0.74 [An (h/z 0 ) - ^ R ] , 


where 

I - 4.7 (h/L) 

- 0.352 (h/L ) 3 - 1.43 (h/L ) 2 
- 2.22 (h/L) 

- 6.35 (h/L) 


for (h/L) 0, 

for -2 < (h/L) < 0; 
for (h/L) > 0, 




0'74/^n - 1 

0.1326 - 2.341 (h/L) 3 1.278 (h/L) 
- 0.2879 (h/L) 


for — 0.08 ( (h/L) £ 0 , 


for - 2 < (h/L) < - 0.08 


In the formulas for ^ and ip R , approximate curve fits have been substituted 
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for the unstable cases (h/L < 0). [The third-order approximations for * was 

111 

devised by Schultz (1979).] 

The Monin-Obukhov length is determined from the local Richardson number 
(at z = h) as 


z 

L 


R i /0.95 
R t /(1 - 5 R ± ) 
38.227 - 463.71 
308.49 - 3323.9 


R ± + 1442.2 rJ 
R £ + 9010.6 R^ 


for R ± < 0 
for 0 < R a 
for 0.1674 
for 0.1875 


< 0.1674, 

< R t < 0.1875, 

< R t < 0.2. 


The first two approximations are from Haltiner and Williams (1980), the latter 
two are curve fits determined from the relationship 


” ( Z /L) 


The eddy diffusion for temperature and moisture substance, ICj,, i s set 
equal to zero at the surface. Hence, the subgrid fluxes of vapor, 
temperature, etc., are not allowed through the ground surface. The effect of 
evaporation from a rain-soaked ground is parameterized by adding a source term 
(at the lowest grid level above the ground) in each the prognostic vapor and 

thermodynamic equations. The formulation for ground evaporation is presented 
in section 4. 

Other boundary conditions at the ground surface are: 


w = 0, 


5Q 


CD 


9z 


0, and 


9z 

9Q 


IC 


9z 


0 , 


0 . 
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The surface boundary condition for pressure is consistent with the vertical 
equation of momentum at z = o; i.e., 

Values at the ground for precipitating moisture substance (rain, snow, and 
hail) are determined by applying upstream time-differencing to the following 
equations : 

aQ RN rr &Q RN 

3t = R 9z * 

9Q SN 5Q SN 

3t S 3z 

!5a. 

at H 3z » 

hence, precipitating hydrometeors fall through the surface boundary at a rate 
determined by their mean terminal velocity. 

Open lateral boundaries 

Computational constraints dictate a requirement for lateral boundaries 
even though no physical counterpart exist. Open lateral boundaries should 
allow the mean flow and superimposed wave modes to pass freely and 
unobstructed. They should also be formulated such that they guarantee 
conservation; i.e., they should not artifically create mass, water vapor, 
etc. The proper formulation of the lateral boundaries is essential for 
accurate simulations. 

An increasingly popular boundary condition for nonperiodic but open flow 
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boundary conditions is the Sommerfeld radiation condition: 


at + c a r 


o. 


( 16 ) 


where $ is any prognostic variable, r is the space coordinate perpendicular 
to the boundary, and C is phase velocity normal to the boundary. Orlanski 
(1976) determined C locally from (16) at an interior grid point then applied 

it to the radiation boundary condition at the following time step; hence he 
assumed 


where N represents the time level and b-1 represents first interior point 
adjacent to the boundary point. Thus, he assumed that the phase speed at the 
boundary is equal to phase speed at the adjacent interior point from the 
previous time step. Orlanski applied this procedure in several cases of two- 
dimensional flow which was governed by a prognostic vorticity equation. He 
found this formulation to work quite well; allowing disturbances to propagate 
through the boundary with minimal reflection and distortion of the interior 
solution. 

In nonhydrostatic primitive equation models the procedure used by Klemp 
and Wilhelmson (1978a), Clark (1979), and Tripoli and Cotton (1980) is to 
apply the radiation boundary condition to the velocity component normal to the 
boundary, and to apply one or more of the following conditions to the 
remaining dependent variables along the boundary: 1) use upwind differencing 

if the normal velocity is directed outward, 2) set the variable equal to a 
fixed reference value, and/or 3) specify the normal gradient equal to zero. 
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The phase velocity in (1) is either specified as in Klemp and Wilhelmson, or 
determined locally as in Orlanski (1976). However, these formulations often 
led to runaway circulations and unrealistic trends in the domain- wide mass 
fields (Clark, 1979; Tripoli and Cotton, 1980, 1982). A new procedure for 
applying the radiation boundary condition in primitive equation models has 
been formulated by Proctor (1985a). This procedure includes the "Orlanski 
radiation boundary condition"; that is Eq. (16) with the phase speed 
extrapolated from the interior as in Orlanski (1976); however, it is in a form 
consistent with the Adams-Bashf orth time differencing scheme. The procedure, 
as outlined below, differs from the conventional approach and is essentially 
free of mass trends and run-away circulations. In fact, a periodic adjustment 
to the domain-wide pressure field [as in Klemp and Wilhelmson (1978a)] is 
never needed throughout any of the simulations. 

At each of the four lateral boundaries, the boundary conditions are as 

follows : 

(1) The radiation boundary condition is applied to the pressure deviation 
and the components of velocity that are tangent to the boundary 
[e.g., Eq. (16) is applied to u and w at the north and south 
boundaries] ; 

(2) The vertical velocity is set equal to zero at its boundary point 
whenever the flow normal to the boundary is directed into the domain; 

(3) The horizontal velocity and pressure can relax to their reference 
values at boundary points where the radiation boundary condition is 
applied; i.e. , 
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is applied to u at the north and south boundaries; 


dv 

at 


e(V 

o 


v) 


is applied to v at the east and west boundaries; and 


5p 

at 


e p 


is applied to p at all four boundaries when the phase speed 

determined by Eq. (16) is directed into the domain. [The value for 
e is 10~ 2 s' 1 .] 

(4) The component of velocity normal to the boundary is determined from 
continuity by assuming that the normal gradient of three-dimensional 
divergence vanishes; thus at the north and south boundaries 


(— ) 


- (OH + i * 

( ax az ^ b-l * 

b 


likewise, on the east— west boundaries 


c— ) = 
dx) b 


/bv aw 

\ a - + * — ) + ®. , , 

by bz ' , b-l 

D 


where ® fe _ 1 is the divergence evaluated at the first interior point; 

(5) The boundary values for the remaining variables are determined from 

upstream time differencing if the flow normal to the boundary is 

directed outward; otherwise if the flow is inward, they are set to 

their reference values; i.e., 0 = 9 . Q * 0 . Q,™ = 0. 

O v ^VO* ^CL) * 

QlC * °» Qr = °» QgN * °» Qr “ °* 


25 


In addition to the boundary conditions, a second-order filter is applied 
to the u, v, w, 9, and Q y fields along the three columns of grid points 
adjacent to each lateral boundary. Application of the filter is necessary In 
order to eliminate 2Ax and 4Ax waves. The Orlanski radiation boundary 
condition is suspected of not being able to handle the fast changing phase 
speeds of the numerically-generated high frequency waves; thus making 
filtering next to the boundaries necessary. 

The above formulation for the lateral boundaries allows the outward 
propagation of disturbances with minimal reflection. But, of course, it 
cannot account for the influence and inward propagation of disturbances, that 
(in the real world) may lie outside of the model domain. 

Top boundary 

At the top boundary the vertical velocity is set equal to zero and the 
potential temperature is held fixed to its reference value. These conditions 
are not unappropriate if the top boundary is chosen at a reasonably high 
altitude. However, an artifact of this choice is the reflection of upward 
propagating gravity waves. To reduce wave reflection, a filter and sponge 
(Perkey and Kreitzberg, 1976) are applied within the four rows below the top 
boundary. 

In the application of the sponge the local rate terms for u, v, w, p, and 
9 are multiplied by a weighting coefficient Wg, which is a function of the 
level beneath the top boundary. The weighting coefficients for w and 9 are: 

0.0 for I = b 

0.4 for I - b-1 

0.7 for I = b-2 

0.9 for I = b-3 

where b refers to the grid points for w or 9 at the top boundary level; 


V r > ’ 
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b-1, the first level beneath the boundary; etc. For u, v, and p a porous 
sponge condition is assumed and the coefficients are given as. 

! 0 .7 for I = b-1 

0.9 for I = b-2 . 

The values of u, v, and p at the top boundary are assigned derivative boundary 
conditions as follows: for horizontal velocity a free slip boundary condition 

is assumed; that is. 


for the pressure deviation, its vertical gradient is assumed to vanish at the 
top boundary, i.e., 


3z 


0 . 


Similarly at the top boundary, the subgrid eddy viscosity and water 
substance variables are 


3z 



= 0 , 


and 


3q cd 8 Qic *r 3q sn 3 % 

3z = 3z " JE Tz w = 0 


The mixing ratio for water vapor, on the other hand, is specified to its 
reference value as Q v = Q vo * 
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4. CLOUD MICROPHYSICS 


The cloud hydrometeors are subdivided into five bulk categories: 1) 

cloud droplets, 2) ice crystals, 3) rain, 4) snow, and 5) hail. The 
hydrometeors comprising the cloud droplets and ice crystal categories are 
assumed to be small and nonprecipitating. The ice crystals are assumed to 
have a monodispersive size distribution. The remaining categories represent 
precipitating hydroraeteors and are assumed to have a continuous inverse 
exponential size distribution. All hydrometeors except ice crystals are 
assumed spherical. The parameterization of the cloud microphysics is similar 
to those described by Orville and Kopp (1977), Lin et al. (1983) and Rutledge 
and Hobbs (1983). 

The size distribution for rain is taken as (Marshall and Palmer, 1948) 


“<V - "or <- VV’ 


(17) 


where N(D R ) is the number of raindrops per unit diameter per unit volume, D R 
is the raindrop diameter, A R is the inverse of the slope of the rain 
distribution, and Nq R is the intercept* Similarly the size distribution for 
snow is assumed as (Gunn and Marshall, 1958) 


S <V ■ 8 os exp ( - VV ; 


(18) 


and the distribution for hail or graupel is taken as (Federer and Waldvogel, 
1975) 


N <V ■ n oh exp w- 


(19) 
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The slope factors can be determined from the above distributions as (e.g., Lin 
et al., 1983) 


a r- '^n* 6 ,> 0,Zi ' 


A s ■ «>oW’"fcs 6 s>°' 25 - 


A H ' <p o < V 1tN OH 5 H )0 ' 25 ’ 


( 20 ) 

( 21 ) 

( 22 ) 


where 6 W > 6 S » and are respectively the densities of water, snow, and hail. 

In the tnicrophysical parameterizations the intercept for rain is assumed 
constant and the intercepts for snow and hail are functions of height only. 

The values for the intercepts are derived from measured size distributions. 

For rain, observations during a thunderstorm by Sekhon and Srivastava (1971) 
lead to a value of N QR = 2.5 x 10 7 m 4 . This value is larger than that found 
by Marshall and Palmer (1948) for widespread rain; but, is nearly identical to 
the N qr for hurricane rainfall reported by Lord et al. (1984). The intercept 
for snow is determined from data reported in Houze et al. (1979) as 

N 0S " 5,5 x 1C)6 ex P 1“ 0-088 ( T 0 " T )J m" 4 . 


where T Q is the temperature of the reference environment (degrees Kelvin) and 
T m is the melting-point temperature (273.16 K). The intercept for hail is 
taken as 


( 4 x 10 4 m 4 

( 4 x 10 4 exp[- 0.088 (T - T )] 


for T < T , 
o — M’ 

for T > T . 
o M 
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The hail intercept at altitudes above the melting level is assumed constant; 
the value of 4 x 10* m - * was deduced by Orville and Kopp (1977) from data 
reported in Federer and Waldvogel (1975). The decrease in N qH with 
temperature at altitudes below the melting level is an attempt to crudely 
approximate the decrease in N qH that is due to the more favorable melting of 
the smaller hail sizes. For example, in a two-dimensional simulation of a 
hailstorm by Kopp et al. (1983), the hail intercept was found to be about an 
order of magnitude lower at the ground than above the melting level. 

The hail category is represented by two different types of spherical ice 
particles: moderate density graupel and hail. As in Cotton et al. (1982), 

only one type is allowed at a given grid point and time. Hail is assumed 
present only when the computed mean graupel diameter exceeds 5 mm. The two 
types of particles differ in that the assumed density for graupel particles is 
lower. Only hail particles are allowed wet growth; otherwise, the 
parameterization of the two types of particles are identical. In the 

following text both types of particles will be referred to as hail. 

2 

The density of hail is 

-3 - -3 

( 900 kg m if > 5 x 10 m 

6 = ! _o _ -3 

H ( 450 kg m if D G i 5 x 10 m » 

where the mass weighted mean diameter of the graupel particle is 
D G = 4 [p Q Q h /450 u N oh ]°* 25 [“]• 

2 The value for the density of the hail particles is taken from Vittori and di 
Caporiacco (1959); the assumed value for moderate density graupel is derived 
from data in Pruppacher and Klett (1978). 


30 



The values of the other densities are assumed as 1000 kg m -3 for liquid 

—3 

water and 100 kg m for snow (e.g., Rutledge and Hobbs, 1983). 

Terminal Velocities 


Rain 

The mass-weighted mean terminal velocity for rain can be determined as 
(Rutledge and Hobbs, 1983) 

77 C N( V " <V "r <V % 

W R° 7 . (23) 

S 0 N<D r ) M (D r ) dD R 

where ^(D^) is the terminal velocity of a raindrop having diameter D d and 

K 

3 

mass M(D^) = ^ D R /6. Rutledge and Hobbs' (1983) polynomial fit to the 

experiment data of Gunn and Kunzer (1949) yields (all units MRS) 

Wr(°r) = “ 0.267 + 5150 D R - 1.0225 x 10 6 D R + 7.55 x 10 7 D R . (24a) 

Also from the data of Gunn and Kinzer, an alternative approximation for 
terminal velocity is (Liu and Orville, 1969) 

VV ‘ 843 D r°' 8 ' (2«b) 

Eq. (24b) is less exact than (24a) but leads to a more simple integration of 
the microphysical equations. Substitution of Eqs. (24a) and (17) into (23) 
yields a mass-weighted mean terminal velocity for rain as 
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(25a) 


4 7 2 

W = [- 0.267 + 2.06 x 10 A D - 2.045 x 10 A 
R 1 R K 

+ 9.06 x 10 9 A^ 3 ] (1.2 /P 0 )°' 4 

(for A > 1.3132 x 10~ 5 m). Substitution of Eqs. (24b) and (17) into (23) 

R 

yields 

W R = 843 r(4.8) A r ' 8 /6 (25b) 

Values for the mean mass-weighted terminal velocity are computed from Eq. 

(25a) since it is more precise and computationally efficient. A correction 
factor is included in (25a) in order to account for the change in fallspeed 
with air density (Foote and Toit, 1969). Eqs. (24b) and (25b) are used only 
in the development of the microphysical parameterizations . 

Snow 

The terminal velocity of snow, in nature, varies slowly with increasing 
particle size; more significant variances in fallspeed are due to snow type 
(e.g., rimed dendrites; graupel like snow). The fallspeed assumed in the 
model is deduced from the data of Locatelli and Hobbs (1974) and Jiusto and 
Bosworth (1971). The terminal velocity for snow is assumed to vary only with 
air density, and is given by (where units are in the MRS system) 

W_ - 1.1 (1.2/p )°* 5 . (26) 

Hail 

From McDonald (1960), the terminal velocity of a hailstone having 
diameter is 
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(27) 


VV -I'* 5 , V 3 p„ c d1 0 ' 5 - 

The drag coefficient for hailstones is assumed constant with a value of 
C]) = 0.45 (Macklin and Ludlam, 1961). The mean mass-weighted terminal 
velocity for hail is determined by substituting Eqs. (19) and (27) into an 
expression similar to (23) and integrating over all particle sizes; giving, 

»H* 1 ’ 09375 t 4 ’ 8 V 3 p 0 c ol°* 5 a h°' 5 - < 28 > 

Parameterization of Microphysical Production Terms 

For most of the microphysical interactions, the production terms are 
parameterized by integrating over the assumed size spectrum. For example, if 
the mass growth of a single raindrop due to a particular interaction is dM/dt, 
then the rate of production of Q R is given by 

dQ oo 

PRODUCTION RATE = ^ M N (D R ) ^ 

P o o 

A description of the microphysical production terms which appear in Eqs. 
(5) - (11) follows below. The production terms have units of per second; 
unless otherwise stated all units are expressed in the MRS system. 

Condensation and evaporation of cloud droplets 

Condensation of water vapor into cloud drops is assumed to occur at rate 
which maintains saturation with respect to water vapor. Likewise, evaporation 
of cloud droplets is assumed to occur at a rate such that either saturation is 
maintained or all available cloud droplets are depleted. From Proctor (1982), 
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\ if + Q CD > °’ 

-OCD ” \ + ^CD < ° 

where At is the time step and \ is defined: 

X = (Q - Q )/[ 1 + Q L 2 /C R T 2 ]. 

V ^sv ; 1 sv V p V 

[A derivation of the above equation is in Appendix D.] The saturation mixing 
ratio for vapor^ is determined from 


PCDWV1 = — 


Q 


sv 


e 

sv 


[Q v + ^]/p, 


(30) 


where, e is the ratio of gas constants (e - “ 0.622), and e g ^ is the 

saturation vapor pressure with respect to water. 


Mean cloud droplet radius 

A mean radius for cloud droplets is needed in several microphysical 
parameter izations . It Is determined by assuming that the number of cloud 
droplets per unit volume (n CD ) remains constant; thus, the average mass of a 
cloud droplet is 

S CD = Q CD p o /n CD* (31> 

and the mean radius of a cloud droplet is 


"^The customary definition for saturation mixing ratio is Q gv - e gv £/[p — e sv ] 
(e.g.. Berry et al., 1945). However a reexamination shows that 
q =’p /p = e e/(p-e ) = e [Q + e]/p. In firestorm simulations the 

customary definition was lound to give negative saturation mixing ratios in 
areas of high temperature, where e gv exceeded p. 
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r CD = f 3 Q CD P o /4 71 6 W n CD ] 


(32) 


1/3 


The value assumed for n CD can be estimated from the condensation nuclei 
spectra, or if not available, climatology of the area being modeled. For 
example, typical values in extreme continental areas are ~ 10 9 m -3 ; while in 
contrast, values in maritime regions can be as low as ~ 10^ m -3 . Note that 
for a given concentration of cloud droplet water, Eq. (32) implies that the 
mean cloud droplet radii are larger in maritime clouds than in continental 
clouds. 

Cloud Ice crystals 

The treatment of ice crystals follows that of Rutledge and Hobbs 
(1983). The ice crystals are assumed to be hexagonal plates and to have a 
monodispersive size distribution. The number concentration of ice crystals is 
assumed to be given by the concentration of ice nuclei active at temperature T 
(e.g., Fletcher, 1962): 

9 

10 T < 230.95, 

Pic exp[0.6 (T m - T)] for T M > T >_ 230.95, 

0 T > t m » (33) 

where n IC is the number of ice crystals m~ 3 , and P IC is a constant usually 
* 2 3 

taken as 10 m“ . An arbitrary upper limit of 10 9 crystals nf 3 is assumed. 

No ice crystal processes occur (except for instantaneous melting) at 
temperatures above 273.16 K. 
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Ice crystal initiation 


Following Rutledge and Hobbs (1983), the Ice crystal process is initiated 
by assuming the immediate presence of small ice crystals having an initial 
diameter of 12.9 pm, whenever the air is saturated with respect to ice. The 
rate of production of ice crystals is computed at temperatures less than 
268.16 K and is given by 

PICWV1 = MIN [10* 12 n IC /p o , \ g ] ( 34 ) 


where 


X 

s 


5 (Q V ” Q si )l1 + 


si 


2 2-1 
L/C R T Z ] . 
s p v 


The saturation mixing ratio with respect to ice is 


Q 


si 


'si 


(Q v + 0/ p » 


(35) 


where e g ^ Is the vapor pressure with respect to ice. 

The lesser of the two rates in Eq. (34) is chosen so as to guarantee that 
the computed growth of the ice crystals will not exceed the vapor available 
for growth. The latter rate has been replaced by a more exact formulation 
than that assumed by Rutledge and Hobbs. The new formulation accounts for the 
implicit adjustment of the saturation mixing ratio (see Appendix 0). 


Deposition and sublimation of ice crystals 

The mass rate of change due to the depositional growth of an ice particle 
(e.g., Pruppacher and Klett, 1978) is (where A* and B* are expressed in terras 
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of the reference environment) 


dM/dt - C[Q v - Q^IFj/IQ^A* + B*J; (36 ) 

where , 

* L v L 

A 1 ~ k_T ^r""T and B i = (P n D „) 1 » 

T o v o u w 

where kj is the thermal conductivity of air and D w is the diffusivity of water 
vapor in air. 

For a small hexagonal plate-like ice crystal the capacitance coefficient 
c ia 4 D IC ; the ventilation factor is assumed unity (F J = 1); and the mean 
ice crystal mass is Q I( , P 0 /n IC (e.g., Rutledge and Hobbs, 1983). 

Substitution into Eq. (36) yields 


PI = 4 n 


IC D IC [Q v" Q si ]/p o [Q si A 1 + B il. 


where the mean diameter of the hexagonal plate-like ice crystals is given by 
Rutledge and Hobbs as 


D I(f 16 * 3 Wic p o /n icl 


0.5 


(37) 


If supersaturation with respect to ice exists (i.e., Q > q ) then 

V S1 


PICWV2 = MIN [PI; X /Atl, 

s 


(38a) 


or if subsaturation exists (i.e., Q < Q ) then 

v si 
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PICWV2 = MAX [PL, - Q IC Mt, X g /At] ; 


(38b) 


The smaller of the two arguments are taken as the rate In (38a) so as to 
guarantee that the computed growth of ice crystals by deposition does not 
exceed the vapor available for growth. Similarly in (38b), the greater of the 
three arguments are taken so as to guarantee the computed sublimation rate 
does not exceed the vapor deficit, nor exceed the amount of ice crystal water 
available. 


Growth of ice crystals due to riming 

The fall velocities of ice crystals and cloud droplets are very small and 
can be ignored in many cloud processes. However, the difference in the fall 
speeds between cloud droplets and ice crystals is important when ice crystal 
growth due to the accretion of supercooled cloud droplets is considered (i.e., 
riming). The growth equation due to riming for an ice crystal (Orville and 
Kopp, 1977) is 


2 

dM uD 

dt' P <! CD — (AV) E iccd , 


(38) 


where AV is the difference in terminal velocity of the ice crystals and cloud 
droplets, and E ICC0 is the collection efficiency of ice crystals collecting 
supercooled cloud droplets. From the above equation, the rate of production 
due to riming is 


PICCDl = £ Q cd n r( , (AV) E ICCD 


(39) 
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The collection efficiency is assumed zero if the cloud droplet radius, 
r Cn» is -*- ess than 6 x 10 ^ m, or if the ice crystal diameter D is less than 
150 x 10 ^ m; otherwise, a collection efficiency of 0.5 is assumed. These 
values assumed for E ICCD were derived from the data presented in Pitter 
(1977). His theoretical results showed that collection between cloud droplets 
and ice crystal plates would not occur if the radius of the droplets are less 
than 6 ym, or if the plate diameters are less than 150 ym. 

Melting of ice crystals 

In regions where T > 273.16 K ice crystals are instantaneously melted 
into cloud droplet water, hence the rate of production of ice crystals due to 
melting is 

PICCD2 = - Q T /At. 

Spontaneous freezing of cloud droplets 

Spontaneous freezing of cloud drops occur if T < 233.16 K ; that is 

PICCD3 = Q CD /At. 

Autoconversion of cloud droplets into rain 

In nature, slowly falling cloud droplets often collect other cloud 
droplets; and by this process, may eventually grow into raindrops. The 
parameterization of this process, which results in the conversion of cloud 
droplet water into rainwater, is called the autoconversion of rainwater 
processes. Kessler (1969) hypothesized an autoconversion rate which was 
initiated when the cloud droplet water exceeded a certain threshold value. 
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Other formulations developed by Berry (1968), Simpson and Wiggert (1969), and 
Berry and Reinhardt (1974b) are more elaborate than the Kessler formulation, 
since they are based on data obtained from detailed stochastic growth 
models. The Berry-Reinhardt formulation should be superior to the Berry 
formulation since the former is based on more recent experiments which utilize 
an improved numerical treatment. 

The Berry-Reinhardt formulation is assumed here; it is determined from 
the expression for the average autoconversion rate as (Pruppacher and Klett, 
1978) 


PRCD1 = l 2 / T 2 if L 2 > °» and T 2 > ° ; < 40 > 

where from Berry and Reinhardt (1974a, 1974b), L 2 and T 2 ate 

L 2 = 0.027 [100 r’ 3 r ^ - 0.4] Q^, (41) 

T 2 = 3.72 [r’ - 7.5f 1 (p^f \ < 42 > 


and 


r' = 10 6 (o/ 0.38) 1/3 r^. (43) 

These parameters are based on the evolution of a bimodal droplet spectrum from 
an initially unimodal (gamma) droplet size distribution; specifically, r' is a 
droplet radius parameter, T 2 represents the time required (in seconds) for the 
predominant radius of the larger mode 4 to reach 50 pm, and L 2 is the liquid 

4 The mode which has the larger sized droplets. 
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water content (g g *) at time T 2 associated with the larger mode. 
Autoconversion is not computed unless both L 2 and T 2 exceed their threshold 
values . 

Threshold values of cloud droplet water needed before autoconversion can 
take place, as based on Eqs. (40) - (43), are given in Table 1. According to 
this formulation autoconversion is unlikely in extreme continental conditions; 
rain, if it does occur, must be initiated by the Bergeron process. This is 
supported by observations in continental areas such as the High Plains, which 
typically find that rain is initiated by the melting of snow, graupel, and 
hail (e.g., Dye et al., 1974; Cannon et al., 1974; Knight et al. , 1974; 

Heymsf ield, 1982). On the other hand, in maritime clouds, the autoconversion 
of rain from cloud droplet water is important, if not essential, for rain 
formation. For example, "warm” rain clouds (convective rain clouds whose tops 
never penetrate above the melting level) are often found within maritime 
regions. Autoconversion of rain from cloud droplet water is the only process 
by which rain can be initiated in "warm" clouds. Note from Table 1 that the 
threshold for autoconversion is very small in maritime conditions. Also from 
Table 1, it is interesting to note that the value assumed by Kessler (1969) 
for the threshold of autoconversion (0,5 g wT J ) corresponds to typical values 
of n^ D and a - 


Collection of cloud droplets by rain 

The growth rate of a single raindrop due to the collection of small cloud 
droplets along its path is determined from the continuous collection equation 
(e.g., Liu and Orville, 1969); i.e., 


dM 

dF- «co W -r Vd- 
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TABLE 1. THRESHOLD CLOUD DROPLET WATER CONTENT NEEDED FOR AUTOCONVERSION 

TO RAIN 


Threshold Cloud Droplet 
Water (g m J ) 

[Derived from Eqs. 

(40) - (43)] 

n CD 

[ # / cm 3 ] 

a 

[Dispersion 

Coefficient] 

Location & Reference 

0.1 

50 

0.366 

Maritime - Simpson & 
Wiggert (1969) 

0.5 

200 

0.30 


1.5 - 2.1 

689 - 927 

0.35 - 0.38 

Upwind of St. Louis - 
Fitzgerald & Spyers- 
Duran (1973) 

2.8 - 3.7 

1157 - 1472 

0.30 - 0.32 

Downwind of St. Louis - 
Fitzgerald & Spyers- 
Duran (1973) 

9.2 

2000 

0.146 

Extreme Continental - 
Simpson & Wiggert (1969) 

1.8 - 17.7 

760 - 3166 

0.12 - 0.32 

Colorado High Plains - 
Knight et al. (1982) 


Following Liu and Orville the integration of the above equation over all drop 

sizes with N(D r ) drops as defined in Eq. (17) and W R (D R ) as defined in Eq. 
(24b); gives 


PRCD2 


30-n; 

76 


a r q cd w r n or e rcd * 


(44) 


The mean collection efficiency of raindrops accreting cloud droplets 
(1 RCD ) ls deter ®ined from a least-square curve fit of experimental data that 
is tabulated in Mason (1971). The formulation, which is based on a constant 
raindrop radius of 1000 pm, is given as 


E RCD^ r CD^ = “ °* 275 44 + 0.26249 x 10 6 r^ 

1.8896 x 10 r^p + 4.4626 x 10^^ r^^. (45) 

which ls computed for: 1.2 x lcf 6 < 7^ < 20 x 10‘ 6 . The collection efficiency 

is set equal to unity, if the mean cloud droplet radius as determined by Eq. 
(32) exceeds 20 x 10*« m (20 p m); if , on the other hand> the cloud drQplet 

radius is less than 1.2 x lO' 6 m (1.2 pm), accretion of cloud droplets by 
raindrops is not allowed. 


Evaporation of raindrops 

Evaporation of raindrops is computed whenever the air is subsaturated and 
there is an insufficient amount of cloud droplet water to erase the 

subsaturation. From Mason (1971) the rate of evaporation of a raindrop having 
mass M is 
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2 * °R (q v I Q sv> r R 

L Q L 

.. 'ntT IF 


v sv ( 

k_ T C R T 
T v 


( R T ‘ 15 + °„P 


ventilation factor for raindrops is given by Ranz and Marshall (1952) 


1/3 „ 0.5 
F- 1+0.3 ^ R e » 


where the Reynolds number is defined: 


R e = W VV 


where v is the molecular viscosity of air. 
m 

With Eqs. (17), (24b) and (25b), integration over the drop size spectrum 
gives a production term for the raindrop evaporation as 


T[Q v -Q gv ] A R [1 + °* 3179 S M 1/3 (A R W R /v m )0 ’ 51 


PRWV1 = 2uN 


, ( 46 ) 


Q sv [A 2 /T - A 3 1 + T/ °w 


where 


A 2* S P o L V 2/R v k T ,and A 3 = P o W 


The evaporation rate may not exceed the amount of rain available; i.e., 


PRWV1 = MAX [PRWV1 , - Q R /At] 


44 



Conversion of Ice crystals to snow 


Following Rutledge and Hobbs (1983), the conversion of ice crystals to 
snow is computed whenever the ice crystal mass exceeds 9.4 x 1(T 10 kg. Hence 
from Eq. (37), the transfer of ice crystal water to snow occurs at a rate 
which maintains a maximum average ice crystal diameter of 500 pm; i.e.. 


PSIC1 = (Q I( , - a) /At, (47) 

where the conversion threshold is: a = 9.4 x 10~ 10 n /o 

IC' ”o* 

Collection of cloud droplets by snow 

The production of snow due to the accretion of cloud droplets is 
parameterized in the same manner as the collection of cloud droplets by 
rain. The terminal velocity for snow [Eq. (26)] is assumed to be independent 
of diameter; this leads to an accretion rate given by 


P 0 = 0.5 u N Q W E 

2 OS ^CD S ”S bCD ’ 


(48) 


If the temperature is less than 273.16 K the accreted cloud water (P 2 ) is 
converted into snow: 


PSCD1 = P 2< 


(49) 


If the temperature is greater than 273.16 K the collected cloud water is 
transferred to rain: 


PRCD3 = P 2# 


(50) 
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The collection efficiency of snow for cloud droplets E gCD , 


Is assumed 


unity (e.g., Lin et^. al_. , 1983). 


Collection of rain by snow 

If the temperature is less than 273.16 K, snow grows by the accretion of 
rain. The rate at which snow accretes rainwater is 


P3 - /" ' 1 (f R +I > s ) 2 » R <V *„ | »<V » <V -a®** 


where W D » W is assumed. Integration after substituting Eqs. (17), (18) and 
R s 

(24b) yields: 


P3 


0.5 it N QS Q r W r A s (13.92 A r + 4.8 + A g 


A?]E 


RS 


(51) 


Snow is produced from the accretion of rain only if T _< 273.16 K 


PSR1 * P3 


if T < V 


The collection efficiency of snow for rain E RS , is assumed unity (e.g., 
Lin et al., 1983; Rutledge and Hobbs, 1984). 


Collection of ice crystals by snow 

The production of snow due to the accretion of ice crystals is 
parameterized in the same manner as the collection of cloud droplets by 
snow. Ice crystals which only occur at temperatures less than 273.16 K are 

accreted at a rate given by 
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( 52 ) 


PSIC2 = 0.5 TT N o Awe 

OS V IC il s S SIC’ 


The collection efficiency of snow for ice crystals is assumed: 


E SIC^ ~ exp[0.38(T T M >]. ( 53 ) 

Eq. (53) is based on experimentally determined efficiencies that are presented 
in Pruppacher and Klett (1978). Note that ice crystals are less likely to 
stick to snow particles at lower temperatures. 


Autoconversion of ice crystal water into snow 

Ice crystals may grow into snow particles by processes of deposition of 
vapor and accretion of cloud droplets. The initiation of snow due to these 
processes are parameterized in Eq. (47). Snow may also form due to the 
collision and aggregation of ice crystals. Following Lin et al. (1983), this 
process is parameterized by an autoconversion formula as, 

PSIC3 - 10- 3 (q ic - 10- 3 /p o ) E SIC , (54) 

which is computed whenever p Q exceeds 10~ 3 kg m~ 3 . 

O 1C 

Melting of snow 

The formulation for melting of ice particles is discussed in Wisner et 
al. (1972). The formulation assumes that the heat required for melting is 
supplied by the following processes: ( 1 ) the conduction of heat from the air, 
( 2 ) the transfer of latent heat due to condensation of water vapor on the 
surface of the ice particle, and (3) the heat supplied by the collection of 
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rain and cloud droplets. The parameterization for the melting of snow 
utilizes Eqs. (18), (26), (48), (51) and (55), and follows Lin et al. 

(1983). The melting of snow results in the production of rain, and is only 
computed when T > T^. The rate of melting is 

2xN 2 

PSR2 = - fk T (T-T M ) + L v Vo ( V Q ssv )] A S 

f 

x [0.86 + 0.28 sJj /3 (TiW s A s /v m ) 0 * 5 ] 

C (T-V 

- Ji iL [PSRl + P2] - PSR1; ( 55 > 

L f 

where Q ssv is the saturation vapor mixing ratio (with respect to liquid water) 
at the surface of the snow particle; i.e., 

Q - e (T ) [Q + e]/P. (56) 

SSV SV M V 

The melting rate may not exceed the amount of snow available, i.e., 

PSR2 = MAX [PSR2, - Q glJ /At]. 

The last term in Eq. (55) is based on the requirement that raindrops be 
shed at the same rate that they are collected. The net effect of raindrop 
collection by melting snow particles is to further enhance the melting rate. 


Growth of snow by deposition and sublimation 

The depositional growth of a spherical snow particle is given by Eq. (36) 
where C is equal to 2xD s> and with the ventilation factor for snow is given 
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as (Hall and Pruppacher, 1976) 


F g = 0.86 + 0.28 S* /3 Re 0 ' 5 , 


( 57 ) 


where Re is now W„D„/v . 

S S m 

Multiplying Eq. (36) by (18) and integrating over all particle sizes with 
Wg defined in (26) yields 


. 2 1 + o.2i 


P4 


PJQ S1 A*+ B,] 


* 


(58) 


Similar to the formulation of Eqs. (38a) and (38b), the rate of 
production of snow by deposition is 


PSWV1 = MIN [P4, \ /At] if Q >q * 

s — x si * 

the rate of production due to sublimation is 

PSWV1 - MAX [P4, - Q SN /At, X g /At] if Q v < q^. 

Deposition and sublimation of snow are not computed if melting is taking 
place; i.e., 


PSWV1 = 0 


if PSR2 < 0. 


Condensation and evaporation of wet snow 

Condensation or evaporation of snow is assumed to occur only during 
melting. Evaporation may take place whenever snow particles melt in air that 
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is below liquid water saturation. Formulation of the process assumes: (1) 

melting and accretion produces a liquid surface completely surrounding the 
particle; and (2) the temperature of the liquid surface is in equilibrium with 
the snow particle. Hence, evaporation occurs if the vapor pressure of the air 
is less than the saturation vapor pressure at 273.16 K; otherwise, water vapor 

is condensed onto the wet snow particle. 

The growth rate of a wet particle where the liquid is in thermal 

equilibrium with ice is 

— = 2 u D p (Q - Q ) Vs* (59> 

(Jt WO V ssv b b 

Multiplying by Eq. (18) and integrating over all sizes yields 
P5 - 2 * N os D w (Q v " Q ssv ) A S 

x [0.86 + 0.21 sJ /3 (* w s V V m )0 ’ 51, (60> 

The rate of condensation or evaporation is only computed when melting is 
occurring; i.e., PSR2 < 0. 

T f o <0 then evaporation from the wet snow results in the 

i L x x ssv* 

production of water vapor: 

PRWV2 = MAX [P5, - P6] ; 
where P6 = PRCD3 - PSR2 + P3. 

The above formulation limits the rate of evaporation to the amount of 
liquid water available from melting and accretion. In other words, the rate 
of evaporation cannot exceed the rate which liquid water is produced on the 
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snow particle. 


On the other hand, if condensation of vapor takes place on the snow 

particles, then the condensed water is converted (at the same rate) into rain 
water: 

PRWV2 = P5. 


Note that for a sufficiently moist atmosphere, rain due to melting snow 
may consist partly of water actually melted from a snow particle, and partly 
of water that was condensed onto the snow particle. 

Spontaneous freezing of raindrops 

The parameterization of the spontaneous freezing of supercooled raindrops 
follows Wisner et^_al_. (1972). This process is included in the TASS 
formulation, even though Lin et al. (1983) has found this process to be 
secondary to drop freezing due to ice crystal collection. The formulation is 
based on the probability function developed by Bigg (1953) from laboratory 
experiments : 


PF = 1 - exp[- a^xD^/6) t G p (T)]; 


(61) 


where 


G f (T) = 


0 if T > 269.16, 

exp[a 2 (T M -T)] - 1 if T < 269.16. 


The constants and a 2 are based on experimental data; the values assumed 
here are taken from Wisner et al. as 
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100 m 3 s \ and 


a l = 


0.66 K 


-1 


From Eq. (61), Wisner et al. obtained an equation for the number of drops 
frozen per unit volume, N(D^), as 

_ = mc^N^) D R G F (T)/6. ( 62 ) 

The formulation of Eq. (62) Is based on the assumption that the number of 
drops frozen is small compared to the total number of drops within a unit 

volume . 

The production rate at which raindrops freeze per unit volume can be 
obtained from eq. (62) as 

” ■ r ’ - f <V M <V 

'o O 

Substitution of Eqs. (17) and integrating yields 


P7 = 20n aj_ Q r A 3 G p (T). (63) 

In the application of Eq. (63), a maximum of 25% of the available rainwater is 

allowed to freeze per time step. This arbitrary upper limit should be 

applied, since the development of Eq. (63) is based on the assumption that 

only a small portion of the number of drops per unit volume actually freeze. 

If Q > 10~ 4 g g” 1 , then hail is produced from the freezing of 
R — 

raindrops : 


PHRl = MIN [P7, Q r /4 At] . 


(64a) 
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Raindrops associated with water contents less than 10" 4 g g" 1 are assumed to 
be too small to classify as hail when frozen, and instead are converted to 
snow. Hence, if Q R < 10 4 g g -1 , snow is produced by the freezing of rain: 

PSR3- MIN [P7, Q r /4 At]. 4 


Raindrop freezing due to collection of ice crystals 

The collection of ice crystals by supercooled raindrops results in the 
production of either hail or snow. Raindrop freezing due to this process 
usually dominates over spontaneous drop freezing (Lin et al., 1983). 

Following Lin et al., if Q R exceeds 10 -4 g g" 1 , the collection of ice 
crystals by rain results in the production of hail. Since the collection of 
an ice crystal by a raindrop results in both particles changing to hail, two 
production terras are needed: the transfer of rainwater to hail due the 

collection of ice crystals, PHR2, and the transfer of ice crystal water to 

hail due to the collection of ice crystals by raindrops, PHIC1. If q D i s i eS s 

“A 1 R 

than 10 g g , the collection of ice crystals by rain results in both the 

transfer of ice crystal water to snow and rainwater to snow. 

Assuming the continuous collection equation, the rate of collection of 
ice crystal water by raindrops is 


P8 = 


15x 

38 “ 


a r q ic w r n or e ric' 


(65) 


The rate at which rainwater is transformed to hail or snow due to the 
collection of ice crystals is (assuming the hydrometeors are uniformly 
distributed throughout the grid volume) 
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P9 ■ ? ' e ric T 4 VV 6 W I ”r n IC “'V «.» 

K o o 


which becomes. 


P9 = 6.96 * n IC A R Q R W R E RI( ,. 


Thus if Q r > 10' 4 g g" 1 , then PHIC1 = P8, and PHR2 = P9; if Q R < 10" 4 g g"\ 
then PSIC4 = P8, and PSR4 = P9. 

As in Lin et al. (1983), the collection efficiency of supercooled rain 
for ice crystals ^ R j(]» I s assumed unity. 

Freezing of supercooled raindrops resulting from the c ollection of snow 

If the temperature is less than 273.16 K and the rain water content 
exceeds 10 -4 g g _1 , hail is produced by the raindrop collection of snow 
particles. Again the process requires two production terms. 

The production rate of hail from rain collecting snow can be determined 
from the rate at which snow accretes rainwater. Hence from Eq. (51), 


PHS1 - PSR1 


if T < T m and Q R > 10 


The rate at which rainwater accretes snow is 


PHS2 - I is Ers 5 (V D S > 2 W 6 s Z “s N( V "‘V dD R dD S' 


o o o 


Substituting Eqs. (17), (18) and (24B), and integrating yields: 


38“ n or q sn w r a r [25 ° A s /63 + 20 A R A s /7 + a R ] e RS’ 
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which is only computed when 


T < 273*16 K, and Q > 10 ^ g g 

R 


Production of hall from the riming of snow 

The riming of large snow particles, especially in regions of large cloud- 
droplet concentrations, may lead to the initiation of hail. 

Price (1985) and Price et al. (1986) assume that the conversion of hail 
from the riming of snow occurs when snow contents, exceeding 0.1 g kg -1 , are 
in coexistence with supercooled cloud-droplet concentrations of 1 g kg - * or 
more. If these threshold conditions are met, they assume that the snow mass 
converted to hail, per time step, is equal to the mass of the snow particles 
greater than some diameter, D Q . The derivation of this formulation assumes a 
continuous size distribution of snow particles; thus with Eq. 18, the mass of 
snow (per mass of dry air) of snow particles greater than D Q is 

OO 

M s ■ “os 5 s { D s exp( W dD s- 

o 

Price (1985) and Price et al. (1986) obtain a production rate by integrating 
the above equation and dividing by the time step and density of air: 


PHS3 = 


SN 

6 A? At 


exp(- -jp-) [D^ + 3 A g + 6 D q A^ + 6 


V 




(69a) 


In order to increase computational efficiency, a least-square curve fit 
is formulated from Eq. (69a) as 
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0 


if 


A < 3.5 x 10 

o 


-4 


PHS3 = 



[-2.31 x 10 -3 + 114 A g -6.05 x 10 5 A^ + 8.59 x 10® A 3 ] 

if 3.5 x 10 4 < A g < 7.5 x 10 \ 
[-0.35 + 445 A g + 2.95 x 10 5 A 3 - 1.21 x 10® A 3 ] 

if 7.5 x 10' 4 < Ag _< 2.0 x 10“ 3 , 
0.75 if Ag > 2 x 10 -3 , (69b) 


_ O 

where D q in (69a) is assumed to have a value of 5 x 10 m. PHS3 is only 

-3 -1 -4 -1 

computed if T < T^, >_ 10 g g > and >10 g g 


Collection by hail of ice crystals 

The production rate for hail due to ice crystal collection is determined 
from the continuous collection equation when T is less than 273.16 K. The 
parameterization follows Lin et. al. (1983) and is formulated in a similar 
manner as PRCD3 in Eq. (44). The production rate for hail collecting cloud 
ice is 

PHIC2 -fS N or Q IC Aj5 W H E hic . (70) 

During hail dry growth the collection efficiency for hail collecting ice 
crystals is E HIC = 0.3, which is based on the experimental data of Latham and 
Saunders (1970). In the case of hail wet growth, Ejjjq is set to unity. 


Collection by hail of cloud droplets 

The production rate for hail due to cloud droplet collection is 


56 


l 



determined from the continuous collection equation, and is formulated in a 

similar manner as PRCD2 in Eq. (44); the production rate is (Lin et. al., 
1983), 

PHCD1 - ~ N O A 3 u EL 

7 OH XD H H T1CD’ (71) 

The mean collection efficiency for hail collecting cloud drops, E , is 

HCD 

based on Langmuir's (1948) theoretical efficiency for potential flow; it is 
given by 


e hcd 


[K s /(K s + °* 5 )] 2 * 


(72) 


The Stokes number (e.g., Byers, 1965) is: 


K 

s 



4 W. 


? 

r /9 
CD X 


U D V 


(73) 


where y D is the dynamic viscosity of air. The mean Stokes number is defined 
from K g by substituting the mean mass-weighted terminal velocity [Eq. (28)] 
and the mean mass-weighted hail diameter [D R = 4^] into Eq. (73); the mean 
Stokes number is given as 


K 

s 


6 w W H r CD /9 


A 


H’ 


Collection by hail of rain 

The rate at which rainwater is collected by hail is (Wisner et al., 1972) 


PHR3 


Jl OO OO TT 2 O 

p7 o r o * <VV 'V“rI 6 W (D R )N(D H> d[, R< ID „ E „ E - (74) 
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Following Wisner et al. , integration of Eq. (74) is possible if differences in 
terminal velocities are approximated by 


'V "R 1 ■ 1 W H - W r'> 

where W and W„ are respectively given by Eqs. (25a) and (28). Substitution 
R H 

of Eqs. (17) and (19) into (74) and integrating gives 

PHR3 = 0.5 *N qh I W r -W r I QgAjj [10 A R + * A R /V R + A R ] E RR . (75) 

The efficiency of hail for rain, E RR , is assumed unity (e.g., Lin et al., 
1983). 

Collection by hall of snow 

The production rate for hail due to snow particle collection is derived 
similar to PHS2 in Eq. (68); it is given as 

PHS4 - N oh W h q sN V 16 V 3 + 16 Vs ' 5 + A H> E HS- (76) 

For hail dry growth the collection efficiency of hail for snow, E RS , is 
assumed equal to 0.1 (Lin et al., 1983); for hail wet growth and temperatures 
greater than 273.16 K, 1® assumed unity. 


Hail melting 

The production rate for the melting of hail is developed similar to PSR1 
in Eq. (55)- The melting of hail results in the production of rain, and is 
computed when the temperature is greater than 273.16 K. Following Lin et al. 
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(1983), the rate at which hail melts is 


P10 


" o — T“~ [ k T < T-T M> + L ° TT P (Q ~Q )] A^ 
P Q J- M v w o v ssv /J H 

X [0.94 + 0.381 s J /3 (« H A H / v /‘ 5 i 
r [ PHCD1 + PHR3 ] ; 


(77) 


where in the derivation of Eq. (77) the ventilation factor for hail is assumed 
as (Mason and Thorpe, 1966) 


F = 0.94 + 0.33 sy 3 R 1 ^ 2 
H Me* 


(78) 


and where now the Reynolds number is defined as R = W D /v 

e H H m' 

The rate of melting may not exceed the hail available and is only 
computed if melting occurs: 


P10 = MAX [P10, - Q /At) , 

H 

P10 = MIN [P10, 0] . 


The rate at which hailwater is transferred to rainwater is 

PHR4 =■ P10 - PHCD1-PHR3 (79) 

The last two terms in Eq. (79) are due to the shedding of collected rain and 
cloud droplet water. 
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Hall wet growth 

Hailstones grow by deposition and the accretion of liquid and solid 
hydrometeors. Growth by accretion of cloud droplets and rain drops is usually 
dominant and becomes more so as the hail diameter increases. Dry growth 
occurs when the surface of the hailstone is at subfreezing temperatures and 
the accreted liquid hydrometeors freeze quickly, leaving the surface 
essentially dry. During wet growth, only a portion of the collected liquid 
water freezes. Wet growth occurs because heat transfer to the surrounding air 
is insufficient to dissipate the excess heat that is released from the 
freezing of accreted water. Thus, the surface temperature of the hailstone 
rises to 0°C, and a portion of the accreted water is shed as rain. According 
to Musil (1970) this process is modified when the hailstone also accretes ice 
particles. Some of the excess heat contained by the hailstone can be absorbed 
by the cooler ice particles; hence, the hailstone's capacity to freeze 
accreted water is increased. 

Parameterization of the wet growth process is similar to the formulation 

in Lin et al. (1983). First, a production term representing the maximum 

capacity for growth from the accreted hydrometeors is computed (PWET). If the 

rate of accreted liquid water exceeds PWET (the maximum growth rate), wet 

growth is assumed to occur. The growth rate for hail is then given by PWET; 

and the remaining portion of the accreted liquid water is shed as rain. 

However, if PWET exceeds the rate at which liquid water is accreted, then dry 

growth is assumed; and all of the accreted water is transferred into hail 

water. This test for wet growth is only computed when the air temperature is 

between 253.16 K and 273.16 K. Only hail particles are assumed to have the 

potential for wet growth; dry growth is always assumed for graupel particles 
_ -3 

(i.e*, when D < 5 x 10 m). 

G 
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The growth equation for wet growth was formulated by Musil (1970) and is 
similar to the production equation for melting. Assuming the ventilation 
factor given by Eq. (78) and integrating over all particle sizes yields, 

2 * N 0H 

PWET = - [k (T-T ) + L D o (Q -Q ) 1 A 2 

P Q T M v w F o vv v v ssv^ H 

x [0.94 + 0.381 S^Of^/v ) 0,5 ] 

C t (T-T ) 

+ [1 ' — f 1 (P«IC2/E hi + PHS4/E HS ); (80) 

which is only computed if 273.16 > T > 253.16 K, and D > 5 x 10 -3 . 

G 

Wet growth occurs if: 

PWET < PHCD1 + PHR3 ; 
in which case, 

PHR5 = PWET - PHCD1 - PHR3 , 


where PHR5 represents the rate at which accreted water is shed as rain. 
Dry growth occurs if: 


PWET 2 PHCD1 + PHR3; 


in which case 


PHR5 = 0, 
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and all accreted water is consumed by the growing hail. 


Condensation and evaporation of wet hail 

The hailstone is assumed to be covered by a film of water whenever 
melting or wet growth occurs. The temperature of the liquid surface is 
assumed to be in equilibrium with the ice surface; thus, condensation on the 
hailstone occurs if vapor pressure of the air exceeds the saturation vapor 
pressure (for liquid water) at 273.16 K. In which case, water vapor condensed 
on the hailstone is transferred to rainwater. Evaporation, on the other hand, 
occurs if the vapor pressure at 273.16 K exceeds that of the air; in this case 
melt- and shed-water (which would otherwise become rain) are transferred to 

water vapor. 

Similar to Eq. (60) the rate of condensation on wet hail is 

V A 2 

Pll = 2 TT N D (Q -Q ) A u 
rij. i. n « QH u w v^ v ^ssv H 

x [0.94 + 0.381 4 /3( VH /V m )0 ’ 5]; (81) 

which is only computed when the hail or graupel particles are "wet"; i.e., 
when either PHR4 or PHR5 are less than zero; otherwise, Pll is set equal to 

zero. 

Condensation on the hailstones occurs if Q v > Q ssv : the water condensed 

on the hail is tranferred to rainwater; 

PRWV3 = Pll. 

Evaporation from the wet hail takes place if Q v ( ^ssv’ 
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PRWV3 = MAX [Pll, - P12], 


(82) 


where the rate at which liquid water is available for evaporation is 
PI 2 =* - PHR4 - PHR5 


By choosing the maximum of the two rates in Eq. (81), evaporation from the wet 
hail cannot exceed the rate at which excess liquid water is collected and/or 
melted. In other words, evaporation occurs as long as the hailstones are wet. 


Growth of hall by deposition and sublimation 

The production rate of hail due to deposition is formulated similar to 
Eq. (58); and is given by 


PHWV1 


2tc N, 


OH^v^sll A Ht °- 94 + °' 381 ^“hVV"' 3 ! 


>0.5, 


P ot Q si A 1 + V 


(83) 


Deposition (or sublimation) is not computed if the hail surface is wet; i.e.. 


PHWV1 =■ 0 


if PRWV3 - PHR4 - PHR5 < 0. 


Evaporation from wet ground 

The rate of evaporation from the ground is computed as a source term in 
the water vapor and thermodynamic equations. The formulation is a crude 
approximation and assumes the ground surface, when it is wet, to be covered by 
a layer of liquid water. 

The flux of vapor to the air from the wet ground is given by 
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dP 

. dM v 

dt w dz * 


where A is the surface area covered by water, M is the mass of water involved 
in the exchange, K„ is the eddy-diffusion coefficient for water vapor, and 

p is vapor density, 
v 

If p Q Ax Ay A z is substituted for M, and the gradient is approximated 
as dP /dz - [P v (h) - P gv (h)]/h, then the production rate of Q v due to surface 


evaporation is 


PWVG = 


dt 


= s 


<V [Q v (h) " Q sv (h)1 


(84) 


2 h 


where 8 is an empirical constant and h - Az/2. In the absence of 

observations, the value of 8 is assumed unity. 

The rate of ground evaporation is computed in the grid cells adjacent to 
the ground, and is only calculated if the ground is assumed wet. For a domain 
that is assumed stationary with respect to ground, the accumulated rainfall 
must locally exceed one millimeter before ground evaporation is computed. If 
the domain is moving (e.g., translating with a convective storm), then PWVG is 
computed locally, only when the rainfall rate exceeds 25 millimeters per hour. 


Source Terms for the Thermodynamic and 
Mnisfnrp Substance Equation 


The source terms for various microphysical interaction are formally 
listed below. 5 


5 Actual computation is not as simple as summing all of the source terms for 
each prognostic equation. Details of the computational procedure are given in 

Section 6. 
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Thermodynamic equation 


The source terms for the rate of diabatic heating which are used in Eq 
(5) are: 


s v = PCDWV1 + PRWV1 + PRWV2 + PRWV3 - PWVG 

S = PICCD1 + PICCD2 + PICCD3 + PSR1 + PSR2 + PSR3 + PSR4 + PHCD1 + PHR1 
+ PHR3 + PHR4 + PHR5 

S g “ PICWV1 + PICWV2 + PSWV1 + PHWV1. 

Moisture substance equatio ns 

The source terms in Eqs. (6) - (11) are as follows: 
for water vapor 


S 

vap 


PCDWVl - PICWV1 - PICWV2 - PRWV1 - PRWV2 - PRWV3 - PSWV1 
- PHWV1 + PWVG; 


for cloud droplet water 


S CD = PCDWV1 ~ PICCD1 - PICCD2 - PICCD3 - PRCD1 - PRCD2 - PRCD3 

- PSCD1 - PHCD1 ; 


for ice crystal water 


S JC = PICWV1 + PICWV2 + PICCD1 + PICCD2 + PICCD3 - PSIC1 - PSIC2 
- PSIC3 - PSIC4 - PHIC1 - PHIC2; 


for rain water 
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S = PRWV1 + PRWV2 + PRWV3 + PRCD1 + PRCD2 + PRCD3 - PSR1 - PSR2 

R - PSR3 - PSR4 - PHR1 - PHR2 - PHR3 - PHR4 - PHR5; 

for snow water 

S =* PSWV1 + PSCDl + PSIC1 + PSIC2 + PSIC3 + PSIC4 + PSR1 

S + PSR2 + PSR3 + PSR4 - PHS1 - PHS2 - PHS3 - PHS4; 

and for hail water 

S =• PHWV1 + PHCD1 + PHIC1 + PHIC2 + PHR1 + PHR2 + PHR3 + PHR4 + PHR5 
H + PHS1 + PHS2 + PHS3. 

Diagnostic Calculation of Radar Reflectivity 

Radar reflectivity fields can be diagnosed from the simulated moisture 
substance fields, since specific size distributions have been assumed. The 
simulated radar-reflectivity fields are useful in the analysis of model 
results; and most important, they can be used to compare and validate the 

model results with actual radar observations. 

The model radar reflectivity fields are diagnosed from the rain, snow, 
and hail fields. The cloud droplet and ice crystal fields, which are composed 
of relatively small particles, contribute very little to the radar 
reflectivity and can be neglected in most applications. 

For rain the radar reflectivity factor is determined by assuming Rayleigh 

scattering: 

Z R - " N < D R> °R <“V <85) 

o 

Integration after substitution of Eq. (17) into (85), gives 
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( 86 ) 


Z R- '- 2 ^ ““ »o« < 


r 6 ~ 3 , 

[mm m ] . 


For dry snow, the equivalent radar reflectivity factor must be assumed; 


i • e • , 


SDRY 


l K l| 2 5 S " 6 

N( V D s 


w 


substitution of Eq. (18) and integrating gives. 


,2 „ 2 


SDRY 


7.2 x 10 


20 * 1 1 




I 2 6 2 

w 


N os A s 


r 6 -3, 

[mm m ] 


(87) 


For wavelengths employed in weather radars, and for temperatures typical of 
meteorological problems, the dielectric factor for water, Jk^; 2 , is 0.93, and 
the dielectric factor for ice, | ^ | 2 , ls 0.21 (Rogers, 1976). The specific 
density of snow is incorporated into Eq. (87) in order to adjust the 
particle diameter to its melted diameter (e.g., Battan, 1973). 

For wet snow the radar reflectivity factor is simply 


snow 


SWET 


7 - 2 - lo2 ° N os A s 


r 6 -3, 

[mm m ]. (88) 


In the case of dry hail, the radar reflectivity factor is calculated in 
the same manner as dry snow; i.e.. 


J HDRY 


7.2 x 10 


20 


IK, 


,2 P 2 


l*W 


,2 p 2 


H N .7 

n oh a h 


r 6 -3, 

[mm m ] . 


(89) 


An empirical formulation for the radar reflectivity of wet hail which 
includes the effects of Mie scattering has been determined by Smith et al. 
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(1975) as 


Z HtfET 


[7.2 


10 


20 


n oh a h^ 


0.95 



-3 

m 


]• 


(90) 


The final radar reflectivity is computed from the sum of the radar 
reflectivity factors of rain, snow, and hail. The wet snow and hail radar 
reflectivity factors are assumed up to the 0° C level. The wet hail radar 
reflectivity factor is also assumed in the updrafts up to the -8° C level, if 

_ -3 

D > 5 x 10 m. 

G 


Evaluation of Microphysical Constants 

Several constants and variables need to be evaluated in order to solve 

for the cloud microphysics. Many of the physical constants vary only slowly 

with temperature and pressure; and thus, are defined in terms of variables 

from the reference environment (see Appendix C). For example, physical 

constants such as latent heat for vaporization of water (L v ), latent heat for 

fusion of water (L f ), latent heat of sublimation of water (L g ), specific heat 

of water (C w ), specific heat of ice, (C I ), dynamic viscosity of air (p Q ), and 

thermal conductivity (1^) are evaluated as functions of the temperature of the 

reference atmosphere (T q ). Other physical constants such as the molecular 

viscosity of air (v ), the diffusivity of water vapor in air (D w ), and the 

m 

Schmidt number (S M ) are defined in terms of temperature and pressure of the 
reference atmosphere- 

Variables such as saturation vapor pressure are more sensitive to small 
changes and must be defined in terms of local variables. The expression for 
deducing the local values of the saturation vapor pressure with respect to 
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liquid water is obtained by integrating the Claus ius-Clapeyr on equation (e.g., 
Pruppacher and Klett, 1978) and expanding the exponent as 


e 

sv 


(T) = 


e svo (T o ) [ 1 + x+x 2 /2+x 3 /3]; 
x - l v [t-t q ]/r v t o 2 . 


(91) 


The saturation vapor pressure with respect to liquid water for the reference 
environment ( e gvo ) is accurately determined from an empirical relation 
(Appendix D). Similarly, the saturation vapor pressure with respect to ice is 

e si (T) - WV [ 1 + x+x 2 /2 + x 3 /3]; 

X - L s [T-T 0 ]/R v T o 2 . (, 2) 

The solutions of the saturation vapor pressures from Eqs. (91) and (92) are 
accurate for most meteorological problems. The computations of the local 

values for the saturation vapor pressures are fast since no exponentials are 
involved. 

The local temperature, T, which is needed in Eqs. (91) and (92) is not a 
working variable and must be diagnosed. Temperature may be evaluated from 
pressure and potential temperature with Poisson's equation; i.e.. 


T 


9 (P/p 

oo 


(93) 


where P qq - 101325 pa and < r R/C p . A computationally efficient formulation 

which avoids exponentiation can be obtained by a perturbation expansion of Eq 
(93) as 
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The above approximation has reasonable accuracy for most meteorological 


problems . 



5. INITIAL AND REFERENCE CONDITIONS 


The reference environment is assumed unsaturated and steady state; its 
values are derived from a vertical profile sounding — representing the 
hydrostatic environment to be modeled. Either an actual rawinsonde sounding, 
a composite of observed soundings, or a profile predicted by a regional 

hydrostatic model 6 may be used. All reference values are a function of height 
only. 

The values for U Q , V Q , T Q , and Q yo are determined from the input 

sounding, and are defined at each vertical level in the model by using a 

spline interpolation. Once these values are determined, then the remaining 

reference variables are computed; the reference pressure, P Q , is obtained by 

integrating the hydrostatic equation; the reference density, p is solved 

o 

from the equation of state; the reference potential temperature, 9 , is 
determined from Poisson's equation; and the remaining thermodynamic variables 
are diagnosed using appropriate formulas (see Appendix E). 

Basic Initial Field 


The basic initial field is assumed to be horizontally homogeneous and is 
defined directly from the reference values; i.e., u(x,y,z,t=o) =■ U Q (z), 
v(x,y,z,t=o) - V Q (z), w(x,y,z,t=o) = 0, p(x,y,z,t-o) - 0, (x,y,z,fc»o) = Q (z) 
and Q v (x,y,z,t=o) = Q vo (z). The remaining moisture substance fields (i.e., 
^CD’ Qic* Qr* ^s* an< ^ are specified as zero. 

Environmental soundings predicted by the Mesoscale Atmospheric Simulation 
System (Kaplan et al., 1982) may be used to initialize TASS. An advantage of 
a regional-scale model sounding is that it can be generated for almost any 
location and time. 
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The external forcing terms in Eqs. (2), (5) and (6) are defined such that 
the basic Initial field remains steady-state. Thus, the fields cannot depart 
from their Initial values unless a perturbation is added. The forcing terms 

are defined as: 


(5u/5t)* = [5x 13 (t=o)/5z]/p o ; 

(5v/5t)* = [5i; 23 (t=*o)/5z] /p 0 ; 

(59/at)* - [5S(9 o )/5z]/p o ; and 
(5Q v /5t)* = [5S(Q v )/5z]/p o . 

Initial Perturbation Field 


(95) 


(96) 


(97) 


(98) 


Numerical cloud modelers have devised various techniques in order to 
trigger convection in their simulations. The most common technique is to 
apply a moist or dry thermal perturbation. For example, Klemp and Wilhelmson 
(1978a, 1978b) have assumed a 21.6 km diameter initial thermal impulse. 
Although the scale is large and difficult to justify, its application did 
promote the development of a realistic appearing storm within a reasonable 
period of time. Other approaches for triggering the development of convective 
storms are: a superimposed velocity and temperature impulse (Schlesinger , 

1984a), a meso-gamma scale vortex (Proctor, 1983), a heating function (Miller, 
1978), a cooling function (Tripoli and Cotton, 1982), a mesoscale forcing 
function for velocity (Schlesinger, 1984b), random heating function (Hill, 
1974; Yau and Michaud, 1982), and topographic uplift due to an isolated 
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mountain (Liu and Orville, 1969; Clark, 1979). Unfortunately all of the 
approaches suffer from some arbritrariness and a lack of understanding of how 
the mesoscale acts to force convection. Smolarkiewicz and Clark (1985) have 
simulated a cumulus field, by including surface energy and moisture balance 
equations, as well as nonhomogeneous terrain into their model. The simulated 
clouds were initiated by the flow over irregular terrain and the nonuniform 
ground heating. This approach for cumulus initiation is progressing in the 
right direction; but it is not yet practical in many applications. Also, 
mesoscale forcing may be more important than boundary layer forcing for 
certain types of storm development « 

In the TASS model convection can be initiated by superimposing a velocity 
impulse and (or) a thermal impulse onto the basic initial field. The 
formulation of the velocity impulse is modified from Schlesinger (1984a). It 
assumes a cylindrical updraft of radius R y and depth H w , and is consistent 
with the anelastic equation for mass continuity; it is given by 


u'(t=o) 


v'(t=o) 


w'(t=o) 


where 


- j— (x-x ) G dF/dz 
p o 

~ 27T (y_y o ) ® dF/dz 
M o 

f 

M , 2 ^ 9 

— exp(- r /R ) [1 - (r/R ) Z ] if r < R 

, Iq w w w 

0 if r > R 

w 


G (r) = 


*7 9 

exp (- r /R 1 
w 

2 *9 

R w exp 


if r < R 

— w 

if r > R , 

w* 


(99) 


( 100 ) 


( 101 ) 
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F(z) - 6.75 z (H w - z) 2 /H^, 


r (x,y) 



+ (y - y c ) 


2 


and 


M = 


P W 
o max 


The maximum updraft speed, W max , occurs at x Q , y 0 , and z - H w /3« 

A thermal perturbation field may be added to the base state as 

0 

T = AT [1 _ (r/R 0 ) 2 - (2z-Z 0 ) 2 /Z 2 ]; (102) 

o 

Q '(t=o) = MAX [0, T ] • 

The thermal perturbation field is an ellipsoid with a maximum horizontal 
radius R 0 and a depth Z 0 . The maximum temperature perturbation, AT, occurs 

at x o» y o and 2 = Z Q^ 2 ’ 
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6. NUMERICAL PROCEDURE 


Choice of Finite-Difference Approximations 


Factors which influence the choice of finite-difference approximations 
are accuracy, economy, and long-term stability. Conservative schemes should 
be used, otherwise artificial generation of mass, momentum, vorticity, and 
energy would obviously invalidate the solutions. Especially desirable are 
schemes for space derivatives introduced by Arakawa (1966), which obey certain 
integral constraints on quadratic quanities, such as kinetic energy. These 
schemes are termed quadratic conservative or energy conservative, since when 
applied to advection they conserve both the first and second statistical 
moments of the dependent variable. These schemes are reasonably efficient and 
are especially popular in long-term integrations since they retard if not 
eliminate the development of nonlinear instability 7 . These numerical schemes, 
however, only possess their quadratic-conservative properties in the absence 
of time-differencing errors. An economical time differencing scheme which is 
complimentary to the quadratic-conservative space differencing is the second- 
order Adams-Bashforth method. Lilly (1965) has shown that the Adams-Bashf orth 
method does not artifically generate kinetic energy when used with quadratic- 
conservative schemes; in addition, it has comparable accuracy, yet is more 
efficient when compared to certain second-order iterative methods. The Adams- 
Bashforth method also has comparable accuracy to the Leapfrog time- 
differencing scheme without the problems of time splitting instability; also 


^ t0 the flnite differ encing of nonlinear terms may lead to 
catostrophic rises in variances associated with the shortest resolvable 

wavelengths (e.g., Haltiner and Williams, 1980.) This so-called nonlinear 
instability cannot be eliminated by reducing the time step. 
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the method is not numerically diffusive, as in the case of the first-order 
upstream differencing scheme. Deardorff (1973), has found the Adams-Bashf orth 
method to be more preferable in three-dimensional turbulent boundary layer 

simulations than the popular Leapfrog scheme. 

A significant reduction in the run time of a compressible formulation can 
be achieved with the time-splitting integration procedure. This scheme has 
been developed for cloud models by Klemp and Wilhelrason (1978a) and Cotton and 
Tripoli (1978). The time-splitting procedure results in a substantial savings 
in computer time, yet results in little loss of accuracy, when compared to 
ordinary methods of compressible integration. In this scheme the higher 
frequency terms given by the LHS of Eqs. (2) and (4) are integrated with a 
time step compatible with the propagation of acoustic modes. The remaining 
terms in Eqs. (2) and (4) along with (5) - (11) are integrated with longer 
time steps which are appropriate for anelastic or incompressible flow. 

In the TASS model the time-splitting integration procedure is used, and 
local time derivatives are approximated by the second-order Adams-Bashforth 
method. Space derivatives are approximated by second-order central 
differences in quadratic-conservative form. Details of the numerical 
formulation are given in the following sections. 


Grid 


The variables are arranged on a conventional staggered grid, often 
referred to as the Arakawa C mesh (Haltiner and Williams, 1980). All 
variables other than velocity are computed at a common point within the center 
of each grid cell. At the midpoints of the faces of the grid cells, the 
velocity component normal to the faces is computed. The grid arrangement 
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allows the use of quadratic conservative schemes and has improved accuracy 
over most other grid arrangements (e.g., Haltiner and Williams, 1980). 

Vertical Stretching 


A vertically stretched grid is obtained by continuously mapping the 

actual vertical coordinate z into the stretched vertical coordinate z'. The 

equations used for the transformation are the same as those used by Wilhelmson 
and Chen (1982). 

The vertical coordinate z is mapped into z* as 


z = ( c i + C 2 z ')z ». 


(103) 


A constant grid interval Az' in z' space is determined from Eq. (103) as 

Az " <-* + <[X 2 + * T l/C 2 )°' 5 l/QtL + 2); (104) 

where X = C l /2C 2- Z T is the height of the domain, and NL is the number of 
levels above the ground. The actual height of each grid point can be 
determined from Eq. (103), where: 


z • = (1-2) A z » 
z • = (1-3/2) Az- 


for w at I = l, 2,...,NL + 2 

for all other variables at I = 1,2,...,NL + 2. 


The mapping factor G is determined as 


G = dz'/dz = [C 1 + 2 C 2 z'J _1 ; 


(105) 
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the vertical derivatives can be transformed into z space as 


d/dz = G d/dz' , and d 2 /dz 2 = G (G -£r) . 

When grid stretching is applied the values typically assumed for constants C L 
and C 2 are 0.168 and 6.4 x 10" 6 m, respectively. These values are taken from 
Wilhelmson and Chen except that the value for C-, is an order of magnitude 
greater. This larger value for C 2 results in a more modest stretching; 
resulting in approximately a factor of 5 increase in vertical grid size from 
bottom to top. No stretching occurs (i.e., z - z ) with 1 an 2 

Severe stretching (a large increase in grid size from top to bottom) is not 
recommended with the current turbulence closure scheme. 

The vertical stretching of the grid mesh gives increased vertical 
resolution near the ground at the expense of resolution near the top of the 
domain. A primary reason for including vertical stretching in a cloud model 
is so that downdraft outflows and accompanying low-level features can be more 
adequately simulated. 


Finite Difference Equations 


Time derivatives 

A generalized form of the Adams-Bashforth time differencing, which allows 
for a variable time step (Ochs, 1975) is employed in all large time-step 
calculations as 

Q 1 * 1 = Q N + N (Q) 
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where the operator <j>( Q) i s expressed as 


*"(q) - At B [(i + ) <fa, 

N-l 


N 

► 

L 


At„ N-l 

N-l (9t L 


N 
2At 


], 


(106) 


where the subscript L refers to derivatives taken over the large time step 
The time levels are defined according to the following notation: 


Q N = Q(t), 

Q^ 1 = Q(t + At N ), 

and 

Q N_1 = «<t - Atg.j). 


The u component of velocity at small time level ml is approximated as 


n+1 


u n + f , /9u" _ ,9u" 1 N , 

2m 3t^ ^ 3t^ ^ + m ^ 


(107) 


where there are m small time steps per large time step. Note that if 
N = „n+£ . N+1 n+£+m 

, then u - u .The subscript s signifies that the derivatives 
are taken over the small time step. 

Both the v and w components of velocity, as well as the pressure 
deviation, are approximated in a similar fashion as u in Eq. (107). 


Space derivatives 

The finite-differencing for the space derivatives used in this study are 
expressed in the operator notation of Shuman (1962) and Lilly (1964) as 

6 X Q [Q(x + 7p) - Q(x - |£)]/Ax 
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and 


Q X 5 [Q(x + 4^) + Q(x - ^)]/2 , 


where x may represent y, z, z' or x. 

In the notations described above the nonacoustical terms In Eq. 

(4) can be expressed as 


N 


^ )« - 6 (u 

at 'l x 


u x U X ) - 6 y (u y v X ) - G6 z ,(u Z ’w X ) 


+ u$ X + fV Xy - 


xz' 
v xy. 


5u x * 


+ 2 6 x (K m 6 x u) + 6 y [K M ( 6 y u + 6 x v)1 
+ G6 z , [K* XZ ’ (G6 z ,u + 6 x w)]/p q - ^(K^) - (^r) 

+ v X [6 u + 6 u + G6 f (G6 ,u f )], 

1 w W Z 


* N 

p- ) T = - 6 (v 

at 'l v 


t 


7 X u y ) - 6 y (v y v y ) - G6 z ,(v w y ) 


+ y - fu’ xy 




+ 6 x [K M Xy( 6 y u + 6 x v)1 + 2 6 y ( Vy v) 

+ G^ Z .[K* XZ ( 6 y w + ^ z . v >]/Po ” T 5 y^ K M $ ^ ~ 

+ v y [&_v + 6 yy v + G6 z ,(G6 z ,v')], 


XX 


= - 6 (W x u a ) - 6 <w 

^ ' v Y 


N 


at 


— * 2 * — Z * 

V W ^ 


- y - z 1 v 

t.t J ai 1 “■ 


+ w® Z '+ f^*” + g(H Z - 1 ) 

+ 6 x [K m XZ ’ (G6 z' U + 6 x W)] + V^ yZ ’ ( V + ^' V)] 

OP 

+ 2G6 z ,(K* G6 z ,w)/p o - — 6 z ,(K*p) 

+ v Z f6 w + 6 w + G6 ,(G6 ,w)], 

1 xx yy z z 


(2) and 


(108) 


(109) 


( 110 ) 
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(Ill) 


at 


N 

\ “ " 6 x ( P X u) " 6 y(P y v) - G6 Z ,(P Z w) 

— z * 

+ P$ - P 0 gw 


+ 6 


* (K M 6 x p) + 6 y< K M y « y P) + 0S ,.l^‘ , « I P]/p o , 


where 


$-6u+6v+G6 ,w 

x y z f 

f = 2 Q sin 
f j~ 2 Q cos <t» A , 

and 

* 

K = p K. 
o 

Any of the moisture substance equations can be expressed as 
dO N 

at \ = “ 6 X (Q X u) " 6 y(Q y v) - G6 z ,(Q Z ' W ) 

+ - G6 z ,(Q 2 W Q p o )/p o 

+ 6 X ( V 6 x Q) + VV 6 y Q) + G 6 *' < ^f Z,G6 B Q)/ Po + s ’ 

+ V[6 XX Q+ V +G6 Z ’ (G6 Z* Q, >1 (112) 

where W Q is the terminal velocity of any precipitating moisture substance 
variable Q. Obviously, the equations for water vapor, cloud droplets, and ice 
crystals will not have a term for the terminal velocity. The finite- 
difference equation for potential temperature has a form similar to Eq. (112). 

To guarantee linear numerical stability for advection, an additional term 
has been added to Eqs. (108) - (110), and (112). In this term the numerical 
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diffusion coefficient, v, is defined from criteria presented in the following 
section . 

The local rate of deformation squared which is used in Eq. (12), is 
expressed as 


| DEF | 2 = 2[(6 v u) 2 + (6 y v) 2 + (6 z w) 2 ] 


+ (6 u Xy + 6 v Xy ) + (G6 z ,v yZ ' + 6 w yZ ) 


, 2 


+ (6 w XZ ’ + G6 , u XZ ' ) 2 ~\$ 2 - 


The acoustically-active terms in Eq. (2) and (4) are computed at each 
small time step and are expressed in finite differences as 


* n 
ou, 

at ; 

s 


5 v 

at 


) = 


9w n 
dt s' 


* n 

®£) = 
5t ; 

s 


- H X (6 X P)/P 0 > 

- H y (6 y P)/P 0 » 

- H Z '(G6 z ,p)/p o , 

- tl P «» 


where $ is derived from the current values of u, v, and w; and P is 
approximated from values at the latest large time step, i.e., p 


Orlanskl boundary condition 

The Orlanski radiation boundary condition is expressed in a form 
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consistent with the Adams-Bashforth nethod and Is applied to v and p at 
each small time step as 


(Q 


n+1 


- V 


3C 


,n-l 


At 


wr'-O' 


where At g = At^/m; and Ar is the grid increment of the coordinate normal to 
the boundary. The phase velocity can be obtained from Eq. (16) as 


C n _1 ( Q n-2 ~"- 2 


b-1 - Q^: 1 1 )/At o 

3 <1 - O** 


Which is limited by conditions on the south and west boundaries as 


0 

c n if 

- 0.5 Ar/At 


C n > 0 

0 C n > ~ 0.5 Ar/At 

s 

C n < - 0.5 Ar/At , 
s’ 


and on the north and east boundaries as 


C" > 0.5 Ar/At 

s 

0 < C n < 0.5 Ar/At 

8 

C n < 0 . 

The Orlanski boundary condition for w need only to be applied at the large 

time steps. The procedure is similar to the above, with the large time step 
and time levels substituted instead. 


0.5 Ar/At 


„n 


C 

0 


if 


83 


Numerical Stability Criteria 


Both the small and large time steps and numerical viscosity coefficient 
are chosen so as to enforce the linear stability of the numerical system. The 
critical short time step At* guarantees linear stability for the acoustic 
modes. Linear stability of the lower frequency modes are enforced by the 
critical large time step At* and the numerical viscosity coefficient v. The 
linear stability analysis for the numerical scheme is found in Appendix F. 

The critical large and small time steps are, respectively, 


At£ - 0 . 5 /MAX [ | u| /Ax + | v| /Ay + | W| /Az] , 

-2 “°‘ 5 

At* - 0.5 |tiR MAX(T ) (Ax' 2 + Ay 2 + MAX[Az ])} 
s * ° 

where |w| - MX[|w|, Iw-Sgl. lar « e tl ” e step ls dl '' lded by 

★ 

/2 whenever it is exceeded by the critical time step At^; i.e., 


At 


N4-1 


At N // 2 


if 


At f«-1 > At L- 


The integer number nt small time steps per large time step must be recompute., 
whenever the large time step is changed; it is given by 


m * Integer + 0 ; 


hence the small time step is 


At 

s 


‘W"’ 
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Which is always less than or equal to the critical small time step. 

The numerical viscosity coefficient is determined such that a critical 
amount of diffusion exists in order to counter the small amplification of the 
Adams-Bashf orth method; it is given as 


v (x,y,z, t) 


At 3 [ 


iu.) 

Ax 


+ M + i»J 


Ay 


Az 


/8[Ax 2 + 


Ay 


-2 


Az~ 2 ] 


In cloud simulations v is typically two to three orders of magnitude less 
than the numerical viscosity term is retained in the model formulation, 
but can probably be neglected with no significant impact in most integrations. 
The maximum critical eddy viscosity is defined 

K MAX = t^At^CAx + Ay 2 + Az 2 )] ; 


hence 


V HIN IV Vxi- 

and 

S- MI “ IV W- 

Thls condltlon 8 u arantees linear atability of the diffealon proceea. However 
in typical cloud aimulationa K, and K„ rarely, if ever, exceed K„, v . 
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Storm Tracking 


The numerical domain may be specified as being stationary with respect to 
the ground, or it may translate with a convective system. The latter offers 
advantages in that: 1) truncation error is reduced since the translation of 

fields relative to the domain is minimized; and 2) convective systems, 
especially fast moving storms, can be simulated with a detailed mesh for long 

periods of time* 

The translation speed of the grid is variable and is computed as 
follows. The grid translation components U G and V G are computed every 

x seconds according to: 


U G (t) “ U G 


. . x(t) - X(t-T) j. n „ 

U„( t r) + + U - Z3 


x(t) - x 


_ — y(t) - y 

v„( t > - v_(t-T) ♦ -y^- i> + 0.25 — ~-2 


where 


_ _ fff Qx dx dy dz 

x ” fff Q dx d y dz 


— _ /// Qy dx dy dz 
, and y ~ jy/ Q dx dy dz 


By defining 


2 2 

Q - A Pq w 


where 


A = 


6 + 2.og 


10 


1.0 


C > io 5 s 1 
x, < io “ 5 s -1 , 


and where G = 3v/3x - 9u/3y; the grid will track the cyclonic rotating storm, 

attempting to keep it centered at location (x Q , y Q ) • Typically 
equal to 300 s and (x Q , y Q ) is set to the center location of the grid. 
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Computation of Microphysics 


Procedure 

Computation of the cloud microphysics entails several steps per large 
time step: 

- ■ e P - Tenative values for 0, Q v , Q CD , Q IC , q r , q sn> 
and Q r are computed from Eqs. (5) - (11) in the absence 
of source terms. 


Step 2. T and Q are diagnosed from 0 , p, and e 

svo 

using Eqs. (30), (91) and (94). 


Q v » Qcd» T > and Q sv are adjusted for 


condensation of water vapor and the evaporation of cloud 

o 

drops . 


5 te P 4 - The production of hail from the riming of snow 
is computed and Q H and Q gN are adjusted accordingly. 

S_te£_5. Production of hail or snow due to the collision 
of raindrops with ice crystals and the spontaneous 
freezing of raindrops are computed (see Appendix G); Q R , 
^sv’ and Qsn are sd justed accordingly. 


— 1 _^P 6 ‘ Q sv i is computed from Eqs. (35) and (92). 


The saturation mixing ratios Q and Q c . are adjusted 
in Appendix D. 


using formulas derived 
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Step 7. Initiation of ice crystals, melting of ice 

crystals, freezing of cloud droplets, growth of ice 
crystals due to deposition and riming, and the conversion 
of ice crystals to snow are computed; Q v , Q 

g 

T, Q sv » and Q sv ^ are adjusted accordingly. 

Step 8. The terminal velocities for rain and hail and 
the remaining microphysical interactions are computed. 

Step 9 . The potential temperature is adjusted for the 
latent heat released. 


Numerical seeding 

"Numerical seeding" in simulated clouds may occur due to the spurious 
presence of very small values of Q^, Qgfl» or Qh* ot ^ er words, very small 

values of rain, snow, or hail, artifically produced by truncation error or 
boundary reflection, may grow in areas where it should not occur. 

In order to prevent numerical seeding the following procedure is used 
following step 1 in the microphysical procedure: 


Q 

Q 

Q 


CD = Q CD 

+ V 

and 

1r- 0 

if 

IC = Q IC 

+ ^SN* 

and 

Q s»- 0 

if 

IC = Q IC 

+ q h , 

and 

Q h - o 

if 


o < q r < io -7 , 
o < q s(1 < io ' 7 

0 < Q H < 10 -7 . 


In other words, rain is converted to cloud droplets if the rain water content 
is positive and not less than 10 ^ g g ; and etc. 

This procedure also increases the computational efficiency of the model, 
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since some of the microphysical interactions need not be computed at grid 
points where small values of Q R , Q gN , or Q H formerly existed. 

Negative water 

In nature Q y , Q^, Q IC , Q R , Q gfJ , and are positive definite quantities; 

but because of numerical approximations, may be negative at some locations in 

the computational domain. In this model's numerical framework, relatively 

small amounts of negative water are produced, and are treated periodically as 

described below. At the end of every 50 time steps the following procedure is 

applied to each of the moisture substance variables; i.e., 0=0 0 0 

» ^ x v’ "CD * "ic* 

Qr» ^SN* and ^H : 

1) The total integrated mass of Q in the domain is 
computed; i.e., 

TMQ = /// Q dx dy dz; 

2) next, negative values at each grid point are set 
equal to zero: 

Q - 0 if Q < 0; 

3) the total integrated mass of Q after eliminating the 
negative values is computed: 

PTMQ * /// p Q dx dy dz; 
o 
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4) the values of Q are multiplied by a factor, rq so 
that the integrated mass is the same as TMQ; i.e., 


Q - Q r q> 


where J 


r Q = (TMQ + gx w )/(PTMQ + » w > • 

Hence, the procedure eliminates the negative values of Q and approximately 

conserves the total integrated mass of Q. 

Negative values of water substance are treated only periodically since 
retention of negative values may result in more accurate space differencing. 
In computing the microphysical production terms, negative values of Q are 

treated as zero. 


Model Code 

The model code is written in Cyber FORTRAN 200, using 64-bit word 
lengths. With the exception of much of the microphysics, the code is almost 
completely vectorized; vectorization is programmed explicitly, along each 
horizontal plane. Configuration of the model domain is made flexible; the 
grid sizes, the number of vertical levels, and the number of grid cells in 
either the x or y direction are input parameters. The model code has a 
restart capability, and simulated data from selected fields is transferred as 

9 The threshold constant a w is given a value of 1<T 8 times the horizontal area 
of the domain. 
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output to disk at a selected time inverval. 

The program runs in-core on the NASA Langley VPS-32 supercomputer. For a 
63 x 63 grid with 32 vertically-stretched layers, the memory required is 5.5 
million words, and the ratio of simulation time to computation time in cloud 
simulations is roughly 2:3. For a 43 x 43 grid with 27 vertically-stretched 
layers the memory required is 2.4 million words, and the ratio of simulation 
time to computation time is roughly 2:1. Actual run time depends on many 
factors such as grid size, number of grid points, speed of environmental 
winds, as well as the intensity and area of cumulus convection. Computations 
with vertically-stretched grids may be several times slower due to the time 
step constraint by the small mesh size near the ground. 

A timing algorithm was run over several time steps during the mature 
phase of a simulated supercell storm. Table 2 shows that roughly half of the 
computational time is used in computing the microphysics. The actual 
percentages may vary with case, time, and the number of grid points assumed. 
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Table 2. BREAKDOWN OF COMPUTATIONAL TIME 
ASSUMING A 63 x 63 x 33 GRID 


Computation of 


Percent of Time 


Microphysics 

Velocity and Pressure (including boundary conditions) 

Remaining Prognostic Variables (including boundary conditions) 
Subgrid Eddy Viscosity 
Near Boundary Filters 
Other 


52% 

26% 

12 % 

3% 


2 % 


5% 


Total 


100 % 
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7. TEST CASES 


Before attempting to evaluate the model against observed data sets, the 
model can be checked out against simplified atmospheric cases. These test 
cases are useful in uncovering any basic flaws in the model coding and 
formulation. Also, these test cases (and additional simulations with the two 
dimensional axisymmetric version of TASS) are useful in evaluating the 
procedure used for the open lateral boundaries. 


No Shear and Unidirectional Shear of the Environmental Winds 

Two simulations are conducted using a composite sounding, representative 
of the atmosphere near Del City, Oklahoma, on 20 May 1977. However, in the 
first case, the environmental winds are removed; hence, convection is 
simulated in a calm environment. In the second case, only the V component of 
the environmental wind is removed (Fig. 2); thus convection is simulated with 
an environment having unidirectional windshear. In both of the cases the 
model domain is 40 km in both horizontal directions and 20 km In the vertical 
direction. The horizontal grid size is 1000 m and the vertical grid size is 
being defined by 31 vertically-stretched layers. Coriolis force 
is neglected. Convection is triggered by assuming a temperature perturbation 
defined by Eq. (102), with AT = 3°C, z Q = 3000 m, and R Q = 10,000 m. The 
perturbation is centered horizontally within the domain at x - o, y - Q . 


Simulation with no shear 

The simplest test case is to assume an axisymmetric perturbation within 
an environment having no ambient wind. It Is expected that the ensuing 
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Fig. 2. Composite sounding for Del City, Oklahoma on 20 May 1977. The wind 
flags represent u component of winds only. Each full barb represents 

5 m s" 1 . 
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convection should be symmetric and stationary with respect to the vertical 
axis of the initial perturbation. 

Figs. 3 and 4 portray y-z cross sections through the axis of the updraft 
at a time shortly after the maximum updraft speed is attained. Note that w, 
v, p, T and radar reflectivity fields all show symmetry about y = o. 

Simulation with unidirectional shear 

In the case of unidirectional shear one expects to see the development of 
counter-rotating vortices which move obliquely from the direction of the wind 
shear. These diverging circulations are associated with the splitting of the 
updraft into two diverging parts (e.g., Wilhelmson and Klemp, 1978b). This 
case is useful in order to determine if the model is capable of simulating 
storm splitting, and is particularly useful in evaluating the storm tracking 
algorithm. 

Results from this simulation are shown in Figs. 5-8. The horizontal 
coordinates along the grids are relative to the position in which convection 

was initiated (x = o, y = o), and change with time due to the translation of 
the domain. 

Pronounced splitting of the updraft is evident by 60 min (Fig. 5b), and 
by 90 min, the northernmost updraft cell begins to exit through the north 
boundary. From Fig. 6b it is evident that the northernmost cell has 
anticyclonic rotation, while its cyclonic-rotating counterpart remains near 
the center of the grid. Each of the two cells are propagating in opposite 
directions normal to the mean tropospheric wind (which is from the west at 8 m 
s ). Hence, the cyclonic cell is propagating to the right of the mean wind, 
and the anticyclonic cell to the left of the mean wind. The cyclonic cell, 
however remains near the center of the grid, since the domain tracking 
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Fig. 3 
case . 
and c) 
in a) 
dashed 





Y (KM) 

Simulated field distributions in y-z plane for the no wind shear 


The fields are a) v component of velocity, b) w component of velocity, 
pressure deviation from environment. The contour intervals are 3ms 
, m s’ 1 in b) and 0.25 mb in c) . Negative values are contoured with 


line . 
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Fig. 4. Same as Fig. 3 except that the fields are a) temperature deviation 

from the environment and b) radar reflectivity. The contour intervals are 3°C 
in a) and 10 dBZ in b) . 
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Fig. 6. Same as Fig. 5 except that plots represent wind vector fields at a) 
30 min and b) 90 min. The wind vectors are determined from the horizontal 
velocity after the environmental winds have been removed. 
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f simulated clouds for unidirectional 


Fig. 7 . Three-dimensional perspectives o 
wind shear case at a) 30 min, b) 60 min, c) 90 min and 
Perspectives viewed from ESE (vertical coordinate in z 


d) 120 min. 
space ) . 




TIME (MIN) 


PMIN 



B 


TIME (MIN) 


Fxg. 8. Time - z plots of a) maximum pressure deviation and b) minimum 
pressure deviation for unidirectional wind shear case. The contour 
interval is 0.2 mb and values range from 3.8 mb in a) to -3.2 mb in b) . 
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algorithm continually adjusts the translation speed of the domain. After 90 
min, the anticyclonic rotating updraft exists through the north lateral 
boundary with no detrimental consequences to the numerical solution. New 
updraft cells form along the southwestern flank of the old cell (Fig. 5d) as 
rain cooled outflow undercuts the potentially unstable air. Three-dimensional 
perspectives of the simulated clouds (which are determined from the cloud 
droplet and ice crystal fields) are shown in Fig. 7 at 30, 60, 90 and 120 min. 

The minimum and maximum pressure deviation in each horizontal level are 
plotted in Fig. 8 as a function of height and time. No obvious mass trends 
exist in the solution, even when the anticyclonic updraft exists through the 
north boundary after 90 min. The greatest values of maximum pressure 
deviation (< 3.8 mb) occur at the surface after 40 min, and are associated 

with the onset of precipitation. 

This case successfully demonstrates the models' capability of simulating 
storm splitting. Also demonstrated is the ability of the model domain to 
preferentially translate with the convective cell having cyclonic rotation. 

The lateral boundary conditions appear stable and allow the outward 
propagation of convective cells without detrimental effects. 

Lateral Boundary Condition Test 

Several experiments have been performed with both the axisymmetric 
version and 3-D version of TASS in order to evaluate the lateral boundary 
conditions. These experiments have been reported earlier in greater detail in 


Proctor (1985a). 



Two-dimensional axisymmetric simulations 

The lateral boundary condition procedure was tested in the two- 
dimensional model by simulating the outward propagation of a downburst gust 
front. Two model domain sizes were used: a 3 km radius x 5 km deep and a 
5 km x 5 km domain. The grid resolution was 40 m in both the radial and 
vertical direction. The downburst circulation was initiated by specifying a 
distribution of hail along the top boundary (see Proctor, 1985b; and Chuang et 
al., 1984 for additional details). The melting of hail and the evaporation of 
rain cool the air and drive the circulation. A severe test to the lateral 
boundary conditions occurs as the gust front propagates through the lateral 
boundary. 

Results from the two experiments are shown in Figs. 9 and 10. Two 
experiments were conducted with everything identical except the domain size. 
The results from the smaller 3 km x 5 km domain are in the left columns of 
Figs. 9 and 10, while results from the 5 km x 5 km domain are in the right 
columns. Comparison between the two experiments show little difference. The 
greatest error occurs in the pressure deviation field, with maximum pressure 
in the smaller domain being roughly 0.1 mb lower. These experiments 
demonstrate that the lateral boundary condition procedure allows the outward 
propagation of the gust front with only minimal reflections and almost no 
alteration of the interior solution. 

Three-dimensional simulation 

The boundary conditions were also tested in the three-dimensional 
simulation of convection for two different domain sizes. The large domain 
covered a horizontal area of 60 km x 60 km while the smaller domain covered 30 
km x 30 km. Both experiments had equal horizontal grid lengths of 1 km and 
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RADIUS (KH) RADIUS (KM) 


Fig. 9. Field distributions from 2-D axisymmetric model of a) simulated radar 
reflectivity, b) radial velocity, and c) vertical velocity. Plots in left 
column are from the 3 km x 5 km domain simulation, while plots in right column 
are from the 5 km x 5 km domain simulation. The contour interval is 10 dBZ in 
a) and 2 m s ^ in b) and c) . 
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Fig. 10. Same as Fig. 9 except that the fields are d) temperature deviation 
from the environment, e) approximate stream function, and f) pressure 

deviation from the environment. The contour intervals are 1°C in d) and 0.1 
mb in f ) . 
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the domain extended to a depth of 18.75 km with a constant vertical grid 
length of 750 m. No precipitation processes were included in these 
simulations • 

A pair of convective cells were similarly initiated in both experiments, 
and remained near the corners of the 30 km x 30 km domain experiment 
throughout the 30 min simulation. The proximity of the convective cells to 
the corners should provide a good test for the lateral boundary conditions. 

A horizontal cross section of the vector field and cloud boundary at 
z = 8250 m is shown in Fig. 11. The simulations agree very well as evidenced 
in Fig. 11. 


106 




Fig. 11. Horizontal cross section of the cloud boundary and wind vector 
field at z = 8250 m for a) 60 km x 60 km domain, and b) 30 km x 30 km 
domain at t = 30 mins. Winds are relative to the translating grid. 
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8. SUMMARY AND CONCLUSIONS 


In this report the theoretical formulation of the TASS model is 
described. The model contains a nonhydrostatic and compressible equation set 
and ice-phase microphysics, with prognostic equations for momentum, pressure, 
potential temperature, water vapor, cloud droplets, ice crystal, rain, snow, 
and hail. A new lateral boundary condition procedure is used which allows 
minimal distortion of the interior flow and the outward propagation of 
convective cells, without the mass trends which sometimes plague other 
boundary condition procedures. A diagnostic surface boundary layer 
formulation is introduced, which is based on similarity theory. The model 
utilizes the time splitting integration procedure with time derivatives 
approximated by the second-order Adams-Bashf orth method. Space derivatives 
are approximated by second-order quadratic-conservative differences. The 
model includes an algorithm which allows the domain to translate along with a 
convective cell, even at variable speeds. In a storm splitting case it was 
shown that the algorithm allows the domain to translate with the convective 
cell having cyclonic rotation. 

Two philosophical approaches are usually assumed in the parameterization 
of cloud microphysics with bulk models. The approach assumed by the Colorado 
State University modelers (e.g., Tripoli and Cotton, 1980; Cotton et al. , 

1982) is that the slope of the hydrometeor size distributions 

(i.e., y and A^) remain constant and the intercepts (i.e., N qr , N QS , 

and Nq R ) vary with water content. On the other hand, the approach assumed by 
the South Dakota School of Mines (e.g., Orville and Kopp, 1977; Lin et al., 

1983) is that the intercepts remain constant and the slopes vary with water 
content. The parameter izations assumed in the TASS model follow the latter 
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approach, and thus differ from those used by the CSU group. It is felt that 
the approach developed by the North Dakota School of Mines Group is more 
simple and better substantiated by observations. Also included in the 
microphysics of the TASS model is an autoconversion of rainwater formulation 
developed by Berry and Reinhardt (1974b), This formulation, which was 
developed from data produced by a detailed stochastic growth model, has a 
threshold which is dependent upon both the total number of cloud droplets per 
unit volume and the dispersion of the droplet spectrum. In continental areas 
the autoconversion of rainwater may be completely suppressed, as is frequently 
observed. 

The formulation of the TASS model is applicable to a wide range of meso- 
gamma scale and microscale phenomena. The TASS model has been applied to 
experiments ranging from the simulation of simple convective clouds to intense 
firestorms. Model verification, a comparison of TASS simulated results with 
detailed observed data sets, is to be given in a following report. One case 
study has been previously described in Proctor (1985c). In this study, a 
supercell hailstorm which passed through a meteorological data observing 
network is modeled. The quasi-steady structure of the storm was simulated 
throughout 4 V 2 hours of simulation time, and many of the observed features of 
the storm were successfully simulated. 
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Constant in grid stretching function 

Constant in grid stretching function 

Diameter of hexagonal plate-like ice crystal 

Hailstone diameter 

Critical snow diameter in hail conversion process (5 x 10 - ^ m) 
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Snow-particle diameter 

Diffusivity of water vapor in air 

Mean mass-weighted graupel particle diameter 

Mean mass-weighted hailstone diameter 
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Rate of deformation 

Collection efficiency of hail for ice crystals 
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Collection efficiency of hail for snow 
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Collection efficiency of ice crystals for cloud droplets 
Collection efficiency of rain for ice crystals (=1) 

Collection efficiency of rain for snow (=1) 

Collection efficiency of snow for cloud droplets 
Collection efficiency of snow for ice crystals 
Mean collection efficiency of rain for cloud droplets 
Vapor pressure for ice 
Saturation vapor pressure for ice 

Saturation vapor pressure for ice at temperature of reference 
environment 

Saturation vapor pressure for liquid water 

Saturation vapor pressure for liquid water at temperature of 
reference environment 

Saturation vapor pressure for liquid water 
Initial velocity impulse function 
Ventilation factor for hail 
Ventilation factor for ice particles 
Ventilation factor for rain 
Ventilation factor for snow 
Coriolis terms 

Mapping factor for vertical grid stretching 

Function in Bigg’s freezing equation 

Universal wind shear function for heat 

Universal wind shear function for momentum 

_2 

The acceleration due to gravity (9.8 ms ) 

Density ratio term 

Depth of initial velocity impulse 

Height of lowest level of grid points above ground (Az/2) 
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Subgrid eddy viscosity for momentum 

Critical subgrid eddy viscosity for numerical linear stability 
Stokes number 

Subgrid eddy viscosity for heat and moisture substances 

Mean Stokes number 

Dielectric factor for ice (0.21) 

Dielectric factor for water (0.93) 

Mean eddy viscosity within surface layer 
Von Karman constant (0.4) 

Thermal conductivity of air 
Monin-Obukhov length 
Latent heat of fusion for water 
Latent heat of sublimation for water 
Latent heat of vaporization for water 
Parameter in autoconversion formula 
Mass 

Mass of snow per mass of dry air 
Average mass of cloud droplet 

Number of small time steps per large time step 
Time level for large time step 

Number of hail particles per unit diameter Dr per unit volume 
Number of raindrops per unit diameter D R per unit volume 
Number of snow particles per unit diameter D g per unit volume 
Number of vertical grid levels above the ground 
Intercept value in hail size distribution 

Intercept value in raindrop size distribution (2.5 x 10^ m”^) 
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Time level for small time step 
Number concentration of cloud droplets 
Number concentration of ice crystals 
Atmospheric pressure 

Production of cloud droplet water due to condensation 

Production of ice crystal water due to riming 

Melting of ice crystal water into cloud droplet water 

Production of ice crystal water due to freezing of cloud 
droplets 

Rate of ice crystal initiation 

Production of ice crystal water due to deposition 
Probability of spontaneous drop freezing 

Production of hail due to the accretion of cloud droplet water 

Production of hail from ice crystal water due to rain 
collecting ice crystals 

Production of hail due to collection of ice crystal water 

Production of hail due to the spontaneous freezing of raindrops 

Production of hail from rainwater due to rain collecting ice 
crystals 

Production of hail due to accretion of rainwater 
Melting of hail into rainwater 

Rate at which accreted water is shed as rain during hail wet 
growth 

Production of hail due to rain collecting snow 
Production of hail due to rain collecting snow 
Autoconversion of snow into hail due to riming 
Production of hail due to collection of snow 
Production of hail due to deposition 
Autoconversion of cloud droplet water into rain 
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PI 
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Production of rain due to collection of cloud droplets 

Production of rain due to collection of cloud droplets by 

melting snow 

Evaporation of rainwater 

Production of rain due to condensation on melting snow 

Production of rain due to condensation on wet hail 

Production of snow due to accretion of cloud droplet water 

Conversion of ice crystal water into snow 

Production of snow due to collection of ice crystal water 

Autoconversion of ice crystal water into snow 

Production of snow from ice crystal water due to rain 
collecting ice crystals 

Production of snow due to collection of rainwater 
Melting of snow into rainwater 

Production of snow due to the spontaneous freezing of rainwater 

Production of snow from rain due to rain collecting ice 
crystals 

Production of snow due to deposition 

Maximum rate at which hail water can be produced from accreted 
liquid water (used in wet growth calculation) 

Rate at which moisture is evaporated from wet ground 

Depositional growth of ice crystals 

Collection of cloud droplets by snow 

Collection of rain by snow 

Depositional growth of snow 

Condensation on wet snow 

Accretion of liquid water by melting snow 

Spontaneous freezing of raindrops 

Collection of ice crystals by raindrops 
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Collection of raindrops by ice crystals 
Rate of hail melting 
Condensation on wet hail 

Production of liquid water available for shedding 

Pressure of reference environment 

Constant in Poisson's equation (10^ pascals) 

Pressure deviation from reference environment 
Dummy variable 

Mixing ratio of cloud droplet water 
Mixing ratio of hail water 
Mixing ratio of ice crystal water 
Mixing ratio of rain water 
Mixing ratio of snow water 

Sum of mixing ratios of hydrometeor water (e.g., rain, snow, 
hail, cloud droplets, ice crystals) 

The value of Q at its first interior grid point from boundary 

Saturation mixing ratio with respect to ice 

Saturation mixing ratio with respect to ice of reference 

environment 

Saturation mixing ratio with respect to liquid water on surface 
of ice particle 

Saturation mixing ratio with respect to liquid water 

Saturation mixing ratio with respect to liquid water of 

reference environment 

Water vapor mixing ratio 

Water vapor mixing ratio of reference environment 
Dummy variable 

Gas contant for dry air (287.04 J kg * K *) 

Reynolds number 
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Richardson flux number 
Richardson number 

Gas constant for water vapor (461.5 J kg -1 K _1 ) 

Radius of initial velocity impulse 
Radius of initial thermal impulse 
Space coordinate normal to boundary 
Droplet radius parameter 

Radial distance from center of initial impulse 
Subgrid eddy flux of q 

Source term in prognostic equation for cloud droplet water 
Source term in prognostic equation for ice crystal water 

Source term in prognostic equation for hail water 

Schmidt number 

Source term in prognostic equation for rainwater 

Source term in prognostic equation for snow water 

Production of ice from liquid water 
Production of ice from water vapor 
Production of liquid water from water vapor 
Source term in prognostic equation for water vapor 
Atmospheric temperature 
Melting temperature (273.16 K) 

Temperature in reference environment 
Time parameter in autoconversion formula 
Time 

Component of grid translation in x direction 
Westward wind component in reference environment 
Westward wind component 
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Z HDRY 

Z HWET 

Z R 

Z SDRY 

Z SWET 
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= u - U 0 

Horizontal wind component 

Component of grid translation in y direction 
Southward wind component in reference environment 
Southward wind component 

- v " v o 

Weighting coefficient for sponge boundary 
Terminal velocity of hailstones 

Maximum updraft speed of initial velocity impulse 
Terminal velocity of raindrops 
Terminal velocity of snow particles 

Mean mass- weighted terminal velocity for hailstones 

Mean mass-weighted terminal velocity for raindrops 

Maximum magnitude of vertical component of either air velocity 
or mean hydrometeor velocity 

Vertical component of velocity 

West-east coordinate 

Reference point for initial perturbation 
South-north coordinate 

Reference point for initial perturbation 
Radar reflectivity factor for dry hail 
Radar reflectivity factor for wet hail 
Radar reflectivity factor for rain 
Radar reflectivity factor for dry snow 
Radar reflectivity factor for melting snow 
Depth of thermal impulse 
Vertical coordinate 
Aerodynamic roughness height 
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Stretched vertical coordinate 

Constant in subgrid-turbulence formulation (0.35) 

Constant in Bigg's freezing equation (1(T 2 m 3 s" 1 ) 

Constant in Bigg’s freezing equation (0.66 K -1 ) 

17.83786198 

Constant in Fletcher's equation (10 — ^ m~^) 

Subgrid turbulence length scale 

Maximum value of initial temperature impulse 

Difference between fall velocities of cloud droplets and ice 
crystals (0.4 ms ) 

Large time step 

Large time step at level N 

Large time step 

Small time step 

Critical large time step for linear stability 
Critical small time step for linear stability 
Grid increment in x direction 
Grid increment in y direction 
Vertical grid increment 

Constant vertical grid increment in z' space 
Density of hail (900 kg m -3 ) 

Density of snow (100 kg m -3 ) 

Density of water (10 3 kg m -3 ) 

Ratio of gas constants (R/R v =■ 0.62197...) 

Ratio of z/L or vertical component of vorticity 
Ratio of specific heats (C /C =* 1.4) 

p V 7 

Potential temperature 

Potential temperature of reference environment 
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p IC 
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Deviation of potential temperature from environment 

R/C (0.28571...) 

P 

Inverse slope of hailstone size distribution 
Inverse slope of raindrop size distribution 
Amount of vapor available for deposition 
Amount of vapor available for condensation 
Inverse slope of snow particle size distribution 
Dynamic viscosity of air 
Numerical viscosity 
Molecular viscosity of air 
PI (3.1415926. . .) 

Density of (wet) air 

Deviation of air density from the environment 
Partial density of cloud droplet water 
Partial density of ice crystal water 

Partial density of hail water 

Partial density of rain water 

Partial density of snow water 

Density of dry air 

Density of dry air in the environment 
Density of (moist) air in the environment 
Partial density of water vapor 

Partial density of water vapor in the environment 

Partial density of water vapor at saturation with respect to 
liquid water 

Dispersion coefficient for cloud-droplet spectrum 
Turbulence stress or parameter in grid tracking algorithm 
Three-dimensional mass divergence 
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The nondimensional temperature gradient in the surface layer 
The nondimensional wind shear in the surface layer 
Operator in Adams-Bashf orth time difference approximation 
Latitude 

Universal function for heat in surface layer parameterization 
Universal function for momentum in surface layer 
Angular velocity of Earth’s rotation (7.292 x 10 -5 s -1 ) 




APPENDIX B 


Derivation of Equation of State and Equation for Pressure 

Given a volume of air, x, in which the precipitation contained by the 
volume is falling at its terminal velocity, the total mass in x is 


M = M, + M + M , 
d v w 


(B— 1 ) 


where M d is the mass of the dry air in t, ^ is the mass of the water vapor, 
and ^ is the mass of the liquid and frozen water in x. [For precipitation 
falling at its terminal velocity, the drag force exerted by the precipitation 
is in exact balance with its gravitational force; thus the precipitation 
imparts its mass to the air.] 

Dividing (B— 1) by the volume of x gives 


P = P d + p v + p w’ (B-2) 

where p d is the density of dry air, p v is the partial density of water 
vapor, and p w is the partial density of the liquid and frozen water, and p is 
the effective air density; i.e., p - M /x, p - M /x, and p = M /x . Eq . 
(B-2) can be rewritten as 


P - P d [1 + Q v + QJ/RT, or 

p = (P - e v ) [1 + Q v + QJ/RT (B-3) 


B-l 



where Q v = P y /P d > Q w = P w /P d > P is atmospheric pressure, e y is vapor 
pressure, R is the gas constant for dry air, and T is temperature* Since in 


atmospheric problems 


e 

v 



P/e, 


where e is the ratio of gas constants for dry air and water vapor, Eq. (B-3) 
becomes 


p - |t t 1 " °* 608 Q v + V* 


(B-4) 


which is the equation of state modified for the effects of water substance. 
The elastic mass continuity equation expressed in cartesian tensors is 




0 . 


Taking d/dt of (B-4) and substituting into (B-5) gives 


du 

dP ^ PdT_PdG_ p i 

dt * T dt G dt Sx i * 


where G = 1 + 0.608 Q + Q . 

V w 


(B-5) 


(B-6) 


Taking the log-derivative of Poisson’s equation as 


1_ dT _ I d9, R_ i dP 
T dt " 0 dt + C P dt ’ 

P 


where 0 is the potential temperature and C p the gas constant for air at 


B-2 



constant pressure 


Substituting the above equation into (B-6) and rearranging 


terms gives 


I ® 

r) dt 


- P 


dU l P d9 

dx. 9 dt 


P dG 
G dt * 


(B-7) 


where p 


1 - — - C /C . 

Cp V p 


If we assume P(x,y,z,t) * P (z) + p(x,y,z , t), where 5P /d z = - 


p o g > 


then (B-7) can be expressed as, 


5p 

dt 


- u + u 

i dx^ i 


p o 8 6 i3 


+ T) 


P d9 
9 dt 


P dG 
’ 11 G dt 





(B-8) 


where 


far- " p i ^ 


/ (1 + 0.608 Q + Q ), 
v w 


or since Q v and Q w are much less than one, 


pP dG 
G dt 


dQ 

riP [0.608 jjr- + 



] (1 - 0.608 Q v 


V* 
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APPENDIX C 


Constant Stress Layer Approximations for Unsaturated Atmosphere 

The wind shear and temperature gradient in the constant stress layer 
(e.g., Haltiner and Martin, 1980) are respectively 

V 

5V * , 

(c-i) 

T 

90 * . 

5z kz H* (C-2) 

2 2 °* 5 

where V is the horizontal wind speed [V = (u + v ) ] , V* is the frictional 

velocity, T* is the nondimensional temperature, 0 is potential temperature, k 

Is von Harman' s constant, <|> is the nondimensional wind shear, and is the 

nondimensional temperature gradient. Both <t> M and <)»„ are a function of z/L 

M H. 

where L is the Monin-Obukhov length; both <)>,_ and <(>„ are assumed to be 
universal from which the values can be determined from field data. 

Integration of (C-l) and (C-2) from the ground to height h (assuming that the 
surface roughness, z q is much less than h) is 

V - V* G M /k, (C-3) 

and 

9(h) - 0(z=o) = G H /k, (C-4) 

where 
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(C-5) 


«n = *” < h /*o> - *M <r>* 

4 'm e ! W 11 ' ♦*«>!"■ c = h/L! 

o 

G h = 0.74 [in (h/z o ) - 4 > h (£)]i (C"6) 


4> r = f h J\ [1 =* <t> H <C)] dAn C, C = h/L. 
o 


From (C-l) the mean wind shear in the constant stress layer (assuming 


z Q « h) is 


f z V * ^M /kZ dz 

<lr> - —2-^ - v* <yhk. 


dz 


/ k dz 
z o 


(C— 7 ) 


With the aid of (C-3) the above result reduces to 


<~> = V (h)/h. 


( 0 - 8 ) 


Similarly, by integrating (C-2), the mean temperature gradient is 


f z T * V kz dz T G 

<m> = l Zo . 

'dz' f h , kh 

J dz 

z o 


(C-9) 


Since , 


T* = T y V* 2 /g k L, 


where T is the average virtual temperature, and g is the acceleration due to 
v 

gravity; (C-9) with the aid of (C-3) can be expressed as 
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(C-10) 


.59. T v V < h > G H ,h % 

<; az > ,2 _ 2 S/' 

8 h Si 


The stress due to the vertical wind shear in the surface layer is 


% * K M 8? • 


where is the eddy viscosity for momentum. With the aid of (C-l), 
x at z = h can be expressed as 


% ' V h > V * *m<L> /kh - 


Since is independent of height in the surface stress layer, the mean 

stress-layer eddy viscosity, <K M >, can be defined such that 


T z = <V 


thus from (07) and (C-ll), 


<V - % ' <5!» ■ V h > W 




APPENDIX D 


Formulation of the Microphysical Adjustment Scheme 


The development for several formulas which are used in microphysical 
adjustment processes are below. Most of the derivation for the condensation 
adjustment is excerpted from Proctor (1982). 

The methodology for condensation and evaporation follows Asai (1965). 
The condensation adjustment scheme maintains saturation by either the 
condensation of water vapor into cloud droplets, or by the evaporation of 
cloud droplets. However, if during evaporation there is an insufficient 
quanity of cloud-droplet water to maintain saturation, then all of the 
available cloud droplet water is evaporated. The derivation for the 
condensation adjustment is as follows. 

Assuming a pseudoadiabat ic and isobaric process, the saturation mixing 
ratio Q vs » is (e.g. Hess, 1959) 


dQ 


sv 


L Q 
V sv 


dT/T 2 R . 

V 


(D-l) 


Since d9/0 » dT/T in an isobaric process, Eq. (D-l) may also be expressed as 


dQ * L Q d9/R T9. 

SV V SV V 


Now if the saturation surplus is defined as 


(D-2) 


5 ' Q vs = AQ 1 + AQ 2’ 


(D-3) 


D-l 



where, if the air is supersaturated (AQ > 0), the amount of water vapor to be 
condensed is represented by AQ^ . The condensation of AQ^ then releases 
latent heat which increases the air's capacity to store additional water vapor 
by AQ 2 * For the case of cloud droplet water in the presence of subsaturated 
air AQ < 0, the maximum amount of cloud droplet water that can be evaporated 
is AQ^, which is less than AQ since the evaporation of AQ^ reduces the 
air’s temperature and capacity to store water. The warming (cooling) due to 
condensation (evaporation) is then 


A9 ,= L 9 AQ. /C T. 
1 v < 1 p 


(D— 4) 


The increased (decreased) storage AQ 2> as a result of condensational warming 
(evaporative cooling) is obtained by combining (D-2), (D-3) and (D-4), hence 
giving 


AQ 


2 


Q L A9 
sv v 


1 


R T9 
v 


L 2 Q AQ. 
V sv 1 


CRT 
P v 


(D-5) 


Since only a fraction r, of the saturation surplus (deficiency) can be used in 
the condensation (evaporation) process, we may define AQ^ = rAQ. Then by use 
of (D-5), we may solve for r as 


r 


aq l aq x 
AQ~ “ AQ X + AQ 2 


AQ X [AQ X + L 2 AQ 


-1 


. Q /C R T z ] 
1 sv p v J 


or 


[1 + L 2 Q /CRT 2 ] 

1 V SV P V 


-1 


(0-6) 
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Thus the amount of water vapor to be condensed or the maximum amount of cloud 
droplet water that can be evaporated, is given by 


\ = 
v 


rAQ. 


(D-7) 


A similar formulation follows for deposition and sublimation if 
saturation with respect to ice is maintained; i.e.. 


dQ 


si 


L Q .dT/T 2 R , 
s^si v' 


(D-8) 


and 


X » 
s 


( Ui» 


[1 + L 2 Q ./C R T 2 ] 
s ^si p v J 


-1 


(D-9) 


Several useful formulas that are derived from (D-l), and (D-8), which can 
be used to compute adjustments to the saturation mixing ratios are as follows: 


1) the change in Q gv due to ice deposition 

dQ * - L L Q dQ /T 2 R C ; 
sv v s x sv H v' v p 

2) the change in Q gv due to the production of frozen water 
liquid water 


(D-10) 


Qj, from 


3) 


dQ =■> L L,Q dQ T /T 2 R C ; 
sv v f x sv H l' v p’ 

the change in Q gi due to deposition 


(D — 11) 


D-3 



and 4) the change in Q gi due to the production of frozen water Q-p from 


liquid water 



APPENDIX E 


Evaluation of Physical Constants 

Many of the physical constants vary only slowly with temperature and 
pressure. So, in order to simplify the model computations and reduce run 
time, the physical constants are defined in terras of the temperature and 
pressure of the reference environment (T 0 , P Q ), rather than local values. The 
expressions below are very accurate empirical curve fits deduced from 

experimental data. All expressions below are defined in terms of the MRS unit 
system. 

The latent heat for vaporization of water (Pruppacher and Klett, 1978) is 

L ( T ) - 2.50078 x 10 6 (273.16/T ) y , 
v ° o 

where y = 0.167 + 3.67 x 10 - * T . 

o 

The latent heat for fusion of water (Pruppacher and Klett, 1978) for 

T < 273.16 K is 
o — 

L f( T o> = 3-3369 x 10 5 + 2030.6 (T q -273.16) - 10.467 (T Q -273.16) 2 . 

At elevations where T q > 273.16 K, the value for the melting point of the 

latent heat of fusion is assumed; i.e., 

L = 3.3369 x 10 5 f or T > 273.16. 

o 

From the first law of thermodynamics (the conservation of energy) the 

E-l 



latent heat of sublimation for water is 


L (T ) - L (T ) + L (T ). 
so VO f o 

The specific heat of ice (Pruppacher and Klett, 1978) for T q < 273.16 K 


is 


C T (T ) = 2106.1 + 7.327 (T - 273.16) 
I o o 


[Otherwise, for T Q > 273.16 K the specific heat of ice is assumed as: 
C T - 2106.1.] 

The specific heat of water (Pruppacher and Klett, 1978) is 


1 5400 

4217.8 + 0.3471 (T - 273. 16) 2 
° 

4178 + 0.01298 (T - 308. 16) 2 , 

+ 1.591 x 10 (T - 308.16) 
o 


for T < 223.16, 
o 

for 273.16 > T > 223.16, 
— o 

for T > 273.16. 
o 


The dynamic viscosity of air (List, 1971) is 


“d 'V 


1.8325 x 10 


416.16 > T o v 

T + 120 496.16 ' 
o 


The molecular viscosity of air is simply 


v 

m 


(P ,T ) 
o o 




The diffusivity of water vapor in air (Hall and Pruppacher, 1976) is 


E-2 



1.94 


-5, 


D w < P o>V = 2,11 X 10 (101325/P o ) (T q /273.16) 


The Schmidt number is defined as the ratio of molecular viscosity to the 
diffusivity of water vapor; i.e., 


S m < p > T „ ) = V „/ D • 
Moo m w 


The thermal conductivity (Wisner et al., 1972) is 


k T <T 0 ) - 1414 V 


Local values of the saturation vapor pressures are determined from 
functions of the saturation vapor pressure of the reference atmosphere. The 
saturation vapor pressure with respect to liquid water of the reference 
atmosphere (e gvo ), and the saturation vapor pressure with respect to ice of 
the reference atmosphere ( e s j_ Q ) (Wexler, 1976, 1977) are respectively 


e (T ) 
svo o 


= exp[- 2991.2729 T 


-2 


-1 


- 6017.0128 T ~ + 18.87643854 

o 

- 0.028354721 T + 0.17838301 x 10~ 4 T 2 

o o 

- 0.84150417 X 10 -9 T 3 + 0.44412543 x 10 _12 T 4 

o o 

+ 2.858487 ftn T 1 , 
o J ’ 


e , (T ) 
sio o 


= exp[- 5865.3696 T 


-1 


+ 22.241033 + 0.013749042 T 

o 

- 0.34031775 x 10~ 4 T 2 + 0.26967687 x 10 -7 T 3 

o o 

+ 0.6918651 In T 1. 

o J 
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APPENDIX F 


Linear Stability Analysis of Numerical Scheme 
Large time step 

Any of the prognostic equations, excluding acoustical terms, may be 
linearized as 

(F-l) 

The stability analysis 
found if we assume a 



where q = q(x,y,z,t), and u,v,w, and K are constants, 
of the finite difference analogue of Eq. (F-l) can be 
general solution as 


q(x,y,z,t) = Q(t) exp[i(£x + my + kz)] 


(F-2) 


where £, m, and k represent wave numbers, and i is the square-root of minus 
one. 

Following Mesinger and Arakawa (1976), the amplitude factor, X , is 
defined such that 


q N+1 



(F-3) 


where the superscript refers to the time level. 


The solution is 


amplifying 

neutral 

decaying 


if 


m > i 
Ul - 1 
Ixl < 1 



A numerical scheme is considered stable if | A.| _< 1* 

Approximating space derivatives in Eq. (F-l) by second-order central 
differences and the time derivative by the Adams-Bashf or th method; and then 
substituting (F-2) gives, 

\ 2 = X[1 + 1.5 (i5 - |0] + 0.5 (i£ ~ H ), ( F “ 4 ) 


where 


t = At^( | u | /Ax + | v | /Ay + |w[/Az), 
assuming a critical wavelength for advection of 4Ax; and 
p = 4At^K (Ax + Ay 2 + Az 2 ), 

assuming a critical wavelength for diffusion of 2Ax. Solutions for the 
physical mode, and the computational mode, A^* can be obtained from 

(F-4) by using a Taylor-series expansion. The physical mode is 

| X | - 1 + C 4 /4 + l 6 /2 - p - 5^ 2 p/4 + p 2 /2 +. . . .H.O.T. (F-5) 

As p -*■ 0 (advection only) the numerical scheme has a weak instability of 
order £^/4. As £ ■+ 0 (diffusion only) the numerical scheme is conditionally 
stable if p _< 1. 

In order to counter the slight instability of advection, which is maximum 
for 4Ax waves, a numerical diffusion coefficient, v, is defined such that 
| | < 1. Assuming a wavelength of 4Ax for both advection and diffusion, then 
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• » « .H.O.T. , 


1 < |X 1 | = 1 + £ 4 /4 - ji/2 + 


and solving for |i, 




l /2 + ....H.O.T. 


Thus , 


4At L v 


(Ax~ 2 + 


Ay 2 + 


Az“ 2 ) = 


0.5 [At L (|u|/Ax+ |v|/Ay + |w|/Az)]' 


or 


At L [| u| /Ax + 


/Ay + 


w[/Az)] 4 /[8(Ax 2 + Ay 2 + Az 2 )]. (F-6) 


Hence, by defining a numerical diffusion coefficient according to (F-6) and 
setting the following criteria: 


At L = 0.5 [ | u | /Ax + | v | /Ay + | w| /Az] -1 , 


and 


= [4At L (Ax 2 + Ay 2 + Az -2 )]; 


then the linear stability of the numerical system is enforced over the large 
time step. 
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Small time step 


Numerical stability of the acoustically active terms 
linearizing the high frequency terms in Eqs. (2) and (4); 


can be determined by 
i • e • , 


at 


+ a = 0, 

ob Xj 


(F-7) 


and 





bu , 


o 5x, * 
J 


(F-8) 


where P and a = o ^ are both assumed constant- Approximating (F-7) and 
o o o 

(F-8) by the Adams-Bashf or th time difference and central space differences, 
and then climating u gives (where now i, j, and k represent the grid indices 
respectively in the x, y, and z direction) 


rt+l 

P ijk 


n , n-1 

2 P ijk + P ijk 


“ r c [9 


(F-9) 


where 


r = Atqa P /4 = Atr]RT /4, 
c 0 0 o 


and 


ijk 


, n . n 

(p i+ljk + P <- 


i-ljk " 2 P ijk )/Ax + (p i}flk + p i j-lk “ 2 P ijk )/Ay 
+ (p ijk+l + P i jk-1 “ 2 P ijk )/Az * 
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Assuming a general solution as 


p * &(t) exp[i (Ax + my + kz)], 
and substituting into (F-9) gives, 

* + (54 r c - 2)\^ + (1-36 r c )\ + 6^=0. (F-10) 

Eq. (F-10) has 3 solutions which are real. The critical small time step 
is chosen such that | X L | < 1, |\ 2 | < 1, and | \ 3 1 < 1; it is given as 

At g < 0.5 [t)RT q (Ax ^ + Ay ^ + Az ^)] (F-ll) 

Hence, approximation of the high frequency terms in Eqs. (2) and (4) by the 
Adams-Bashforth method and central differences has a conditional numerical 
stability. 
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APPENDIX G 


Raindrop Freezing 

Special treatment is required in computing raindrop freezing, especially 
since the production rates for this process may acquire large magnitudes* 

From Eq • (63) and (64) the rate of depletion of rain due to spontaneous 
freezing is 

DR1A--Q R P7«, . (G-l) 

where P7* = MAX[20 tc A r G p , 0.25/At], and where At is the size of the 
large time step. From Eq. (67) the depletion of rain due to contact freezing 
with ice crystals is 

DRIB =* ~ Q r P9* , (G-2) 

where P9' - 6.96 * n IC A* \ 

Hence, the depletion of rain due to spontaneous freezing and contact freezing 
is 


(dQ R /dt) p = - Q r (P7 ' + P9' ). 

By treating P7' and P9' as constant during the time interval At, the 
depletion of rain during the large time step interval is 


G-l 



(G— 3 ) 


AQr= Q r [1 - exp[- At(P7 ' + P9')]]. 


If p < 0.5, where p = At(P7’ + P9'), Eq. (G-3) can be approximated as 


A Q r = Q r P (1 - 0.5 p). 

The mixing ratios may then be adjusted as follows: 






> 

-4 

-1 

% = q h + 

‘Or 

if 


10 

g g * 






-4 

-1 

^SN = Q SN 

+ % 

if 

Qr 

< 

10 

g g , 


and 

Or 3 * "V 

This procedure prevents the over depletion of rain when the production 
rates for drop freezing are of large magnitude. 
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